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>
23#include <hrleDenseIterator.hpp>
28#include <unordered_map>
32#ifdef LS_TO_VISUALIZATION_DEBUG
33#include <vtkXMLRectilinearGridWriter.h>
38using namespace viennacore;
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;
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;
61 static void removeDuplicatePoints(vtkSmartPointer<vtkPolyData> &polyData,
62 const double tolerance) {
64 vtkSmartPointer<vtkPolyData> newPolyData =
65 vtkSmartPointer<vtkPolyData>::New();
66 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
67 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
68 ptInserter->SetTolerance(tolerance);
70 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
74 polyData->GetBounds(gridBounds);
77 std::vector<vtkIdType> newPointIds;
78 newPointIds.reserve(polyData->GetNumberOfPoints());
79 ptInserter->InitPointInsertion(newPoints, gridBounds);
82 for (vtkIdType pointId = 0; pointId < polyData->GetNumberOfPoints();
84 vtkIdType globalPtId = 0;
85 ptInserter->InsertUniquePoint(polyData->GetPoint(pointId), globalPtId);
86 newPointIds.push_back(globalPtId);
90 newPolyData->SetPoints(newPoints);
93 vtkSmartPointer<vtkCellArray> oldCells = polyData->GetPolys();
94 vtkSmartPointer<vtkCellArray> newCells =
95 vtkSmartPointer<vtkCellArray>::New();
97 vtkSmartPointer<vtkIdList> cellPoints = vtkIdList::New();
98 oldCells->InitTraversal();
99 while (oldCells->GetNextCell(cellPoints)) {
100 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
102 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
105 newCells->InsertNextCell(cellPoints);
108 newPolyData->SetPolys(newCells);
114 newPolyData->GetCellData()->ShallowCopy(polyData->GetCellData());
117 polyData = newPolyData;
122 static void removeDuplicatePoints(vtkSmartPointer<vtkUnstructuredGrid> &ugrid,
123 const double tolerance) {
125 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
126 vtkSmartPointer<vtkUnstructuredGrid>::New();
127 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
128 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
129 ptInserter->SetTolerance(tolerance);
131 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
134 double gridBounds[6];
135 ugrid->GetBounds(gridBounds);
138 std::vector<vtkIdType> newPointIds;
139 newPointIds.reserve(ugrid->GetNumberOfPoints());
140 ptInserter->InitPointInsertion(newPoints, gridBounds);
143 for (vtkIdType pointId = 0; pointId < ugrid->GetNumberOfPoints();
145 vtkIdType globalPtId = 0;
146 ptInserter->InsertUniquePoint(ugrid->GetPoint(pointId), globalPtId);
147 newPointIds.push_back(globalPtId);
151 newGrid->SetPoints(newPoints);
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();
159 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
162 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
170 newGrid->GetCellData()->ShallowCopy(ugrid->GetCellData());
178 removeDegenerateTetras(vtkSmartPointer<vtkUnstructuredGrid> &ugrid) {
179 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
180 vtkSmartPointer<vtkUnstructuredGrid>::New();
183 vtkSmartPointer<vtkIntArray> materialNumberArray =
184 vtkSmartPointer<vtkIntArray>::New();
185 materialNumberArray->SetNumberOfComponents(1);
186 materialNumberArray->SetName(
"Material");
190 vtkDataArray *matArray =
191 ugrid->GetCellData()->GetArray(
"Material", arrayIndex);
192 const int &materialArrayIndex = arrayIndex;
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();
201 for (vtkIdType nextId = pointId + 1;
202 nextId < cellPoints->GetNumberOfIds(); ++nextId) {
204 if (cellPoints->GetId(pointId) == cellPoints->GetId(nextId))
210 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
213 if (materialArrayIndex >= 0)
214 materialNumberArray->InsertNextValue(matArray->GetTuple1(cellId));
219 newGrid->SetPoints(ugrid->GetPoints());
220 newGrid->GetPointData()->ShallowCopy(ugrid->GetPointData());
222 newGrid->GetCellData()->SetScalars(materialNumberArray);
231 template <
int gr
idExtraPo
ints = 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()) {
237 auto &grid = levelSet->getGrid();
238 auto &domain = levelSet->getDomain();
240 int numLayers = levelSet->getLevelSetWidth();
242 vtkSmartPointer<vtkFloatArray>
244 int gridMin = 0, gridMax = 0;
249 for (
unsigned i = 0; i <
D; ++i) {
250 coords[i] = vtkSmartPointer<vtkFloatArray>::New();
252 if (grid.getBoundaryConditions(i) ==
253 Domain<T, D>::BoundaryType::INFINITE_BOUNDARY) {
255 gridMin = std::min(domain.getMinRunBreak(i), infiniteMinimum) -
258 gridMax = std::max(domain.getMaxRunBreak(i), infiniteMaximum) + 1;
262 gridMin = grid.getMinGridPoint(i) - gridExtraPoints;
263 gridMax = grid.getMaxGridPoint(i) + gridExtraPoints;
266 for (
int x = gridMin; x <= gridMax; ++x) {
267 coords[i]->InsertNextValue(x * gridDelta);
273 coords[2] = vtkSmartPointer<vtkFloatArray>::New();
274 coords[2]->InsertNextValue(0);
277 vtkSmartPointer<vtkRectilinearGrid> rgrid =
278 vtkSmartPointer<vtkRectilinearGrid>::New();
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]);
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");
307 viennahrle::ConstDenseIterator<typename Domain<T, D>::DomainType> it(
308 levelSet->getDomain());
311 for (vtkIdType pointId = 0; pointId < numGridPoints; ++pointId) {
315 rgrid->GetPoint(pointId, p);
317 viennahrle::Index<D> indices(grid.globalCoordinates2GlobalIndices(p));
323 if (grid.isOutsideOfDomain(indices)) {
324 indices = grid.globalIndices2LocalIndices(indices);
330 it.goToIndices(indices);
331 if (it.getValue() == Domain<T, D>::POS_VALUE) {
333 }
else if (it.getValue() == Domain<T, D>::NEG_VALUE) {
336 value = it.getValue() + LSOffset;
350 signedDistances->SetValue(pointId, value * gridDelta);
407 rgrid->GetPointData()->SetScalars(signedDistances);
412 void addMetaDataToVTK(vtkDataSet *data)
const {
413 if (metaData.empty()) {
418 vtkSmartPointer<vtkFieldData> fieldData = data->GetFieldData();
419 for (
const auto &meta : metaData) {
420 if (meta.second.empty())
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]);
430 fieldData->AddArray(metaDataArray);
435 WriteVisualizationMesh() =
default;
437 WriteVisualizationMesh(SmartPointer<Domain<T, D>> levelSet) {
438 levelSets.push_back(levelSet);
442 void insertNextLevelSet(SmartPointer<Domain<T, D>> levelSet) {
443 levelSets.push_back(levelSet);
448 void setFileName(std::string passedFileName) {
449 fileName = std::move(passedFileName);
453 void setExtractHullMesh(
bool passedExtractHullMesh) {
454 extractHullMesh = passedExtractHullMesh;
458 void setExtractVolumeMesh(
bool passedExtractVolumeMesh) {
459 extractVolumeMesh = passedExtractVolumeMesh;
462 void setMaterialMap(SmartPointer<MaterialMap> passedMaterialMap) {
463 materialMap = passedMaterialMap;
466 void setWrappingLayerEpsilon(
double epsilon) { LSEpsilon = epsilon; }
469 const std::unordered_map<std::string, std::vector<T>> &passedMetaData) {
470 metaData = passedMetaData;
473 void addMetaData(
const std::string &key,
T value) {
474 metaData[key] = std::vector<T>{value};
477 void addMetaData(
const std::string &key,
const std::vector<T> &values) {
478 metaData[key] = values;
482 const std::unordered_map<std::string, std::vector<T>> &newMetaData) {
483 for (
const auto &pair : newMetaData) {
484 metaData[pair.first] = pair.second;
490 for (
unsigned i = 0; i < levelSets.size(); ++i) {
491 if (levelSets[i]->getLevelSetWidth() < 2) {
492 Logger::getInstance()
494 "WriteVisualizationMesh: Level Set " + std::to_string(i) +
495 " should have a width greater than 1! Conversion might fail!")
500 const double gridDelta = levelSets[0]->getGrid().getGridDelta();
503 std::vector<vtkSmartPointer<vtkUnstructuredGrid>> materialMeshes;
504 std::vector<unsigned> materialIds;
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) {
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));
525 vtkSmartPointer<vtkTableBasedClipDataSet> clipper =
526 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
527 auto topGrid = vtkSmartPointer<vtkRectilinearGrid>::New();
530 LS2RectiLinearGrid(levelSets.back(), 0, totalMinimum, totalMaximum);
536#ifdef LS_TO_VISUALIZATION_DEBUG
538 auto gwriter = vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
539 gwriter->SetFileName(
"./grid_0.vtr");
540 gwriter->SetInputData(topGrid);
542 std::cout <<
"Wrote grid 0" << std::endl;
545 clipper->SetInputData(topGrid);
546 clipper->InsideOutOn();
547 clipper->SetValue(0.0);
548 clipper->GenerateClippedOutputOn();
551#ifdef LS_TO_VISUALIZATION_DEBUG
553 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
554 gwriter->SetFileName(
"./clipped.vtu");
555 gwriter->SetInputData(clipper->GetClippedOutput());
557 std::cout <<
"Wrote clipped" << std::endl;
561 const bool useMaterialMap = materialMap !=
nullptr;
562 materialMeshes.emplace_back(clipper->GetOutput());
563 materialIds.push_back(useMaterialMap ? materialMap->getMaterialId(0) : 0);
565#ifdef LS_TO_VISUALIZATION_DEBUG
567 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
568 gwriter->SetFileName(
"./probed_0.vtu");
569 gwriter->SetInputData(materialMeshes.front());
577 for (
typename LevelSetsType::const_reverse_iterator it =
578 ++levelSets.rbegin();
579 it != levelSets.rend(); ++it) {
580 if (it->get()->getNumberOfPoints() == 0)
585 vtkSmartPointer<vtkRectilinearGrid> rgrid =
586 vtkSmartPointer<vtkRectilinearGrid>::New();
588 rgrid = LS2RectiLinearGrid<1>(*it, -LSEpsilon * counter, totalMinimum,
595#ifdef LS_TO_VISUALIZATION_DEBUG
597 vtkSmartPointer<vtkXMLRectilinearGridWriter> gwriter =
598 vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
599 gwriter->SetFileName(
600 (
"./grid_" + std::to_string(counter) +
".vtr").c_str());
601 gwriter->SetInputData(rgrid);
603 std::cout <<
"Wrote grid " <<
counter << std::endl;
608 vtkSmartPointer<vtkProbeFilter> probeFilter =
609 vtkSmartPointer<vtkProbeFilter>::New();
610 probeFilter->SetInputData(materialMeshes.back());
611 probeFilter->SetSourceData(rgrid);
612 probeFilter->Update();
614#ifdef LS_TO_VISUALIZATION_DEBUG
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());
622 std::cout <<
"Wrote unstructured grid " <<
counter << std::endl;
629 vtkSmartPointer<vtkTableBasedClipDataSet> insideClipper =
630 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
631 insideClipper->SetInputConnection(probeFilter->GetOutputPort());
632 insideClipper->GenerateClippedOutputOn();
633 insideClipper->Update();
635 materialMeshes.back() = insideClipper->GetOutput();
636 materialMeshes.emplace_back(insideClipper->GetClippedOutput());
639 material = materialMap->getMaterialId(counter);
640 materialIds.push_back(material);
645 vtkSmartPointer<vtkAppendFilter> appendFilter =
646 vtkSmartPointer<vtkAppendFilter>::New();
648 vtkSmartPointer<vtkAppendPolyData> hullAppendFilter =
649 vtkSmartPointer<vtkAppendPolyData>::New();
651 for (
unsigned i = 0; i < materialMeshes.size(); ++i) {
654 vtkSmartPointer<vtkIntArray> materialNumberArray =
655 vtkSmartPointer<vtkIntArray>::New();
656 materialNumberArray->SetNumberOfComponents(1);
657 materialNumberArray->SetName(
"Material");
660 materialMeshes[materialMeshes.size() - 1 - i]->GetNumberOfCells();
662 materialNumberArray->InsertNextValue(materialIds[i]);
664 materialMeshes[materialMeshes.size() - 1 - i]->GetCellData()->SetScalars(
665 materialNumberArray);
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);
680 if (extractHullMesh) {
681 vtkSmartPointer<vtkGeometryFilter> geoFilter =
682 vtkSmartPointer<vtkGeometryFilter>::New();
683 geoFilter->SetInputData(materialMeshes[materialMeshes.size() - 1 - i]);
685 hullAppendFilter->AddInputData(geoFilter->GetOutput());
688 appendFilter->AddInputData(materialMeshes[materialMeshes.size() - 1 - i]);
692 auto volumeVTK = vtkSmartPointer<vtkUnstructuredGrid>::New();
693 auto hullVTK = vtkSmartPointer<vtkPolyData>::New();
694 if (extractVolumeMesh) {
695 appendFilter->Update();
699 volumeVTK = appendFilter->GetOutput();
700#ifdef LS_TO_VISUALIZATION_DEBUG
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());
715 removeDuplicatePoints(volumeVTK, 1e-3 * gridDelta);
717#ifdef LS_TO_VISUALIZATION_DEBUG
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);
731 vtkSmartPointer<vtkDataSetTriangleFilter> triangleFilter =
732 vtkSmartPointer<vtkDataSetTriangleFilter>::New();
733 triangleFilter->SetInputData(volumeVTK);
734 triangleFilter->Update();
735 volumeVTK = triangleFilter->GetOutput();
738 removeDegenerateTetras(volumeVTK);
741 addMetaDataToVTK(volumeVTK);
743 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
744 writer->SetFileName((fileName +
"_volume.vtu").c_str());
745 writer->SetInputData(volumeVTK);
750 if (extractHullMesh) {
751 hullAppendFilter->Update();
752 hullVTK = hullAppendFilter->GetOutput();
754 removeDuplicatePoints(hullVTK, 1e-3 * gridDelta);
756 vtkSmartPointer<vtkTriangleFilter> hullTriangleFilter =
757 vtkSmartPointer<vtkTriangleFilter>::New();
758 hullTriangleFilter->SetInputData(hullVTK);
759 hullTriangleFilter->Update();
761 hullVTK = hullTriangleFilter->GetOutput();
764 addMetaDataToVTK(hullVTK);
766 auto writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
767 writer->SetFileName((fileName +
"_hull.vtp").c_str());
768 writer->SetInputData(hullVTK);
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