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>
31#ifdef LS_TO_VISUALIZATION_DEBUG
32#include <vtkXMLRectilinearGridWriter.h>
37using namespace viennacore;
45template <
class T,
int D>
class WriteVisualizationMesh {
46 typedef typename Domain<T, D>::DomainType hrleDomainType;
47 using LevelSetsType = std::vector<SmartPointer<Domain<T, D>>>;
48 LevelSetsType levelSets;
49 SmartPointer<MaterialMap> materialMap =
nullptr;
51 bool extractVolumeMesh =
true;
52 bool extractHullMesh =
false;
53 bool bottomRemoved =
false;
54 double LSEpsilon = 1e-2;
59 static void removeDuplicatePoints(vtkSmartPointer<vtkPolyData> &polyData,
60 const double tolerance) {
62 vtkSmartPointer<vtkPolyData> newPolyData =
63 vtkSmartPointer<vtkPolyData>::New();
64 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
65 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
66 ptInserter->SetTolerance(tolerance);
68 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
72 polyData->GetBounds(gridBounds);
75 std::vector<vtkIdType> newPointIds;
76 newPointIds.reserve(polyData->GetNumberOfPoints());
77 ptInserter->InitPointInsertion(newPoints, gridBounds);
80 for (vtkIdType pointId = 0; pointId < polyData->GetNumberOfPoints();
82 vtkIdType globalPtId = 0;
83 ptInserter->InsertUniquePoint(polyData->GetPoint(pointId), globalPtId);
84 newPointIds.push_back(globalPtId);
88 newPolyData->SetPoints(newPoints);
91 vtkSmartPointer<vtkCellArray> oldCells = polyData->GetPolys();
92 vtkSmartPointer<vtkCellArray> newCells =
93 vtkSmartPointer<vtkCellArray>::New();
95 vtkSmartPointer<vtkIdList> cellPoints = vtkIdList::New();
96 oldCells->InitTraversal();
97 while (oldCells->GetNextCell(cellPoints)) {
98 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
100 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
103 newCells->InsertNextCell(cellPoints);
106 newPolyData->SetPolys(newCells);
112 newPolyData->GetCellData()->ShallowCopy(polyData->GetCellData());
115 polyData = newPolyData;
120 static void removeDuplicatePoints(vtkSmartPointer<vtkUnstructuredGrid> &ugrid,
121 const double tolerance) {
123 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
124 vtkSmartPointer<vtkUnstructuredGrid>::New();
125 vtkSmartPointer<vtkIncrementalOctreePointLocator> ptInserter =
126 vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
127 ptInserter->SetTolerance(tolerance);
129 vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
132 double gridBounds[6];
133 ugrid->GetBounds(gridBounds);
136 std::vector<vtkIdType> newPointIds;
137 newPointIds.reserve(ugrid->GetNumberOfPoints());
138 ptInserter->InitPointInsertion(newPoints, gridBounds);
141 for (vtkIdType pointId = 0; pointId < ugrid->GetNumberOfPoints();
143 vtkIdType globalPtId = 0;
144 ptInserter->InsertUniquePoint(ugrid->GetPoint(pointId), globalPtId);
145 newPointIds.push_back(globalPtId);
149 newGrid->SetPoints(newPoints);
152 for (vtkIdType cellId = 0; cellId < ugrid->GetNumberOfCells(); ++cellId) {
153 vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
154 ugrid->GetCellPoints(cellId, cellPoints);
155 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
157 cellPoints->SetId(pointId, newPointIds[cellPoints->GetId(pointId)]);
160 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
168 newGrid->GetCellData()->ShallowCopy(ugrid->GetCellData());
176 removeDegenerateTetras(vtkSmartPointer<vtkUnstructuredGrid> &ugrid) {
177 vtkSmartPointer<vtkUnstructuredGrid> newGrid =
178 vtkSmartPointer<vtkUnstructuredGrid>::New();
181 vtkSmartPointer<vtkIntArray> materialNumberArray =
182 vtkSmartPointer<vtkIntArray>::New();
183 materialNumberArray->SetNumberOfComponents(1);
184 materialNumberArray->SetName(
"Material");
188 vtkDataArray *matArray =
189 ugrid->GetCellData()->GetArray(
"Material", arrayIndex);
190 const int &materialArrayIndex = arrayIndex;
193 for (vtkIdType cellId = 0; cellId < ugrid->GetNumberOfCells(); ++cellId) {
194 vtkSmartPointer<vtkIdList> cellPoints = vtkSmartPointer<vtkIdList>::New();
195 ugrid->GetCellPoints(cellId, cellPoints);
196 bool isDuplicate =
false;
197 for (vtkIdType pointId = 0; pointId < cellPoints->GetNumberOfIds();
199 for (vtkIdType nextId = pointId + 1;
200 nextId < cellPoints->GetNumberOfIds(); ++nextId) {
202 if (cellPoints->GetId(pointId) == cellPoints->GetId(nextId))
208 newGrid->InsertNextCell(ugrid->GetCell(cellId)->GetCellType(),
211 if (materialArrayIndex >= 0)
212 materialNumberArray->InsertNextValue(matArray->GetTuple1(cellId));
217 newGrid->SetPoints(ugrid->GetPoints());
218 newGrid->GetPointData()->ShallowCopy(ugrid->GetPointData());
220 newGrid->GetCellData()->SetScalars(materialNumberArray);
229 template <
int gr
idExtraPo
ints = 0>
230 vtkSmartPointer<vtkRectilinearGrid>
231 LS2RectiLinearGrid(SmartPointer<Domain<T, D>> levelSet,
const double LSOffset,
232 int infiniteMinimum = std::numeric_limits<int>::max(),
233 int infiniteMaximum = -std::numeric_limits<int>::max()) {
235 auto &grid = levelSet->getGrid();
236 auto &domain = levelSet->getDomain();
238 int numLayers = levelSet->getLevelSetWidth();
240 vtkSmartPointer<vtkFloatArray>
242 int gridMin = 0, gridMax = 0;
247 for (
unsigned i = 0; i <
D; ++i) {
248 coords[i] = vtkSmartPointer<vtkFloatArray>::New();
250 if (grid.getBoundaryConditions(i) ==
251 Domain<T, D>::BoundaryType::INFINITE_BOUNDARY) {
253 gridMin = std::min(domain.getMinRunBreak(i), infiniteMinimum) -
256 gridMax = std::max(domain.getMaxRunBreak(i), infiniteMaximum) + 1;
260 gridMin = grid.getMinGridPoint(i) - gridExtraPoints;
261 gridMax = grid.getMaxGridPoint(i) + gridExtraPoints;
264 for (
int x = gridMin; x <= gridMax; ++x) {
265 coords[i]->InsertNextValue(x * gridDelta);
271 coords[2] = vtkSmartPointer<vtkFloatArray>::New();
272 coords[2]->InsertNextValue(0);
275 vtkSmartPointer<vtkRectilinearGrid> rgrid =
276 vtkSmartPointer<vtkRectilinearGrid>::New();
278 rgrid->SetDimensions(coords[0]->GetNumberOfTuples(),
279 coords[1]->GetNumberOfTuples(),
280 coords[2]->GetNumberOfTuples());
281 rgrid->SetXCoordinates(coords[0]);
282 rgrid->SetYCoordinates(coords[1]);
283 rgrid->SetZCoordinates(coords[2]);
295 auto const numGridPoints = rgrid->GetNumberOfPoints();
296 vtkSmartPointer<vtkFloatArray> signedDistances =
297 vtkSmartPointer<vtkFloatArray>::New();
298 signedDistances->SetNumberOfComponents(1);
299 signedDistances->SetNumberOfTuples(numGridPoints);
300 signedDistances->SetName(
"SignedDistances");
305 viennahrle::ConstDenseIterator<typename Domain<T, D>::DomainType> it(
306 levelSet->getDomain());
309 for (vtkIdType pointId = 0; pointId < numGridPoints; ++pointId) {
313 rgrid->GetPoint(pointId, p);
315 viennahrle::Index<D> indices(grid.globalCoordinates2GlobalIndices(p));
321 if (grid.isOutsideOfDomain(indices)) {
322 indices = grid.globalIndices2LocalIndices(indices);
328 it.goToIndices(indices);
329 if (it.getValue() == Domain<T, D>::POS_VALUE) {
331 }
else if (it.getValue() == Domain<T, D>::NEG_VALUE) {
334 value = it.getValue() + LSOffset;
348 signedDistances->SetValue(pointId, value * gridDelta);
405 rgrid->GetPointData()->SetScalars(signedDistances);
411 WriteVisualizationMesh() =
default;
413 WriteVisualizationMesh(SmartPointer<Domain<T, D>> levelSet) {
414 levelSets.push_back(levelSet);
418 void insertNextLevelSet(SmartPointer<Domain<T, D>> levelSet) {
419 levelSets.push_back(levelSet);
424 void setFileName(std::string passedFileName) {
425 fileName = std::move(passedFileName);
429 void setExtractHullMesh(
bool passedExtractHullMesh) {
430 extractHullMesh = passedExtractHullMesh;
434 void setExtractVolumeMesh(
bool passedExtractVolumeMesh) {
435 extractVolumeMesh = passedExtractVolumeMesh;
438 void setMaterialMap(SmartPointer<MaterialMap> passedMaterialMap) {
439 materialMap = passedMaterialMap;
442 void setWrappingLayerEpsilon(
double epsilon) { LSEpsilon = epsilon; }
446 for (
unsigned i = 0; i < levelSets.size(); ++i) {
447 if (levelSets[i]->getLevelSetWidth() < 2) {
448 Logger::getInstance()
450 "WriteVisualizationMesh: Level Set " + std::to_string(i) +
451 " should have a width greater than 1! Conversion might fail!")
456 const double gridDelta = levelSets[0]->getGrid().getGridDelta();
459 std::vector<vtkSmartPointer<vtkUnstructuredGrid>> materialMeshes;
460 std::vector<unsigned> materialIds;
462 int totalMinimum = std::numeric_limits<int>::max();
463 int totalMaximum = -std::numeric_limits<int>::max();
464 for (
auto &it : levelSets) {
465 if (it->getNumberOfPoints() == 0) {
468 auto &grid = it->getGrid();
469 auto &domain = it->getDomain();
470 for (
unsigned i = 0; i <
D; ++i) {
471 if (grid.getBoundaryConditions(i) ==
472 Domain<T, D>::BoundaryType::INFINITE_BOUNDARY) {
473 totalMinimum = std::min(totalMinimum, domain.getMinRunBreak(i));
474 totalMaximum = std::max(totalMaximum, domain.getMaxRunBreak(i));
481 vtkSmartPointer<vtkTableBasedClipDataSet> clipper =
482 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
483 auto topGrid = vtkSmartPointer<vtkRectilinearGrid>::New();
486 LS2RectiLinearGrid(levelSets.back(), 0, totalMinimum, totalMaximum);
492#ifdef LS_TO_VISUALIZATION_DEBUG
494 auto gwriter = vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
495 gwriter->SetFileName(
"./grid_0.vtr");
496 gwriter->SetInputData(topGrid);
498 std::cout <<
"Wrote grid 0" << std::endl;
501 clipper->SetInputData(topGrid);
502 clipper->InsideOutOn();
503 clipper->SetValue(0.0);
504 clipper->GenerateClippedOutputOn();
507#ifdef LS_TO_VISUALIZATION_DEBUG
509 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
510 gwriter->SetFileName(
"./clipped.vtu");
511 gwriter->SetInputData(clipper->GetClippedOutput());
513 std::cout <<
"Wrote clipped" << std::endl;
517 const bool useMaterialMap = materialMap !=
nullptr;
518 materialMeshes.emplace_back(clipper->GetOutput());
519 materialIds.push_back(useMaterialMap ? materialMap->getMaterialId(0) : 0);
521#ifdef LS_TO_VISUALIZATION_DEBUG
523 auto gwriter = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
524 gwriter->SetFileName(
"./probed_0.vtu");
525 gwriter->SetInputData(materialMeshes.front());
533 for (
typename LevelSetsType::const_reverse_iterator it =
534 ++levelSets.rbegin();
535 it != levelSets.rend(); ++it) {
536 if (it->get()->getNumberOfPoints() == 0)
541 vtkSmartPointer<vtkRectilinearGrid> rgrid =
542 vtkSmartPointer<vtkRectilinearGrid>::New();
544 rgrid = LS2RectiLinearGrid<1>(*it, -LSEpsilon * counter, totalMinimum,
551#ifdef LS_TO_VISUALIZATION_DEBUG
553 vtkSmartPointer<vtkXMLRectilinearGridWriter> gwriter =
554 vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
555 gwriter->SetFileName(
556 (
"./grid_" + std::to_string(counter) +
".vtr").c_str());
557 gwriter->SetInputData(rgrid);
559 std::cout <<
"Wrote grid " <<
counter << std::endl;
564 vtkSmartPointer<vtkProbeFilter> probeFilter =
565 vtkSmartPointer<vtkProbeFilter>::New();
566 probeFilter->SetInputData(materialMeshes.back());
567 probeFilter->SetSourceData(rgrid);
568 probeFilter->Update();
570#ifdef LS_TO_VISUALIZATION_DEBUG
572 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
573 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
574 gwriter->SetFileName(
575 (
"./probed_" + std::to_string(counter) +
".vtu").c_str());
576 gwriter->SetInputData(probeFilter->GetOutput());
578 std::cout <<
"Wrote unstructured grid " <<
counter << std::endl;
585 vtkSmartPointer<vtkTableBasedClipDataSet> insideClipper =
586 vtkSmartPointer<vtkTableBasedClipDataSet>::New();
587 insideClipper->SetInputConnection(probeFilter->GetOutputPort());
588 insideClipper->GenerateClippedOutputOn();
589 insideClipper->Update();
591 materialMeshes.back() = insideClipper->GetOutput();
592 materialMeshes.emplace_back(insideClipper->GetClippedOutput());
595 material = materialMap->getMaterialId(counter);
596 materialIds.push_back(material);
601 vtkSmartPointer<vtkAppendFilter> appendFilter =
602 vtkSmartPointer<vtkAppendFilter>::New();
604 vtkSmartPointer<vtkAppendPolyData> hullAppendFilter =
605 vtkSmartPointer<vtkAppendPolyData>::New();
607 for (
unsigned i = 0; i < materialMeshes.size(); ++i) {
610 vtkSmartPointer<vtkIntArray> materialNumberArray =
611 vtkSmartPointer<vtkIntArray>::New();
612 materialNumberArray->SetNumberOfComponents(1);
613 materialNumberArray->SetName(
"Material");
616 materialMeshes[materialMeshes.size() - 1 - i]->GetNumberOfCells();
618 materialNumberArray->InsertNextValue(materialIds[i]);
620 materialMeshes[materialMeshes.size() - 1 - i]->GetCellData()->SetScalars(
621 materialNumberArray);
627 vtkSmartPointer<vtkPointData> pointData =
628 materialMeshes[materialMeshes.size() - 1 - i]->GetPointData();
629 const int numberOfArrays = pointData->GetNumberOfArrays();
630 for (
int j = 0; j < numberOfArrays; ++j) {
631 pointData->RemoveArray(0);
636 if (extractHullMesh) {
637 vtkSmartPointer<vtkGeometryFilter> geoFilter =
638 vtkSmartPointer<vtkGeometryFilter>::New();
639 geoFilter->SetInputData(materialMeshes[materialMeshes.size() - 1 - i]);
641 hullAppendFilter->AddInputData(geoFilter->GetOutput());
644 appendFilter->AddInputData(materialMeshes[materialMeshes.size() - 1 - i]);
648 auto volumeVTK = vtkSmartPointer<vtkUnstructuredGrid>::New();
649 auto hullVTK = vtkSmartPointer<vtkPolyData>::New();
650 if (extractVolumeMesh) {
651 appendFilter->Update();
655 volumeVTK = appendFilter->GetOutput();
656#ifdef LS_TO_VISUALIZATION_DEBUG
658 std::cout <<
"Before duplicate removal: " << std::endl;
659 std::cout <<
"Points: " << volumeVTK->GetNumberOfPoints() << std::endl;
660 std::cout <<
"Cells: " << volumeVTK->GetNumberOfCells() << std::endl;
661 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
662 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
663 gwriter->SetFileName(
"before_removal.vtu");
664 gwriter->SetInputData(appendFilter->GetOutput());
671 removeDuplicatePoints(volumeVTK, 1e-3 * gridDelta);
673#ifdef LS_TO_VISUALIZATION_DEBUG
675 std::cout <<
"After duplicate removal: " << std::endl;
676 std::cout <<
"Points: " << volumeVTK->GetNumberOfPoints() << std::endl;
677 std::cout <<
"Cells: " << volumeVTK->GetNumberOfCells() << std::endl;
678 vtkSmartPointer<vtkXMLUnstructuredGridWriter> gwriter =
679 vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
680 gwriter->SetFileName(
"after_removal.vtu");
681 gwriter->SetInputData(volumeVTK);
687 vtkSmartPointer<vtkDataSetTriangleFilter> triangleFilter =
688 vtkSmartPointer<vtkDataSetTriangleFilter>::New();
689 triangleFilter->SetInputData(volumeVTK);
690 triangleFilter->Update();
691 volumeVTK = triangleFilter->GetOutput();
694 removeDegenerateTetras(volumeVTK);
696 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
697 writer->SetFileName((fileName +
"_volume.vtu").c_str());
698 writer->SetInputData(volumeVTK);
703 if (extractHullMesh) {
704 hullAppendFilter->Update();
705 hullVTK = hullAppendFilter->GetOutput();
707 removeDuplicatePoints(hullVTK, 1e-3 * gridDelta);
709 vtkSmartPointer<vtkTriangleFilter> hullTriangleFilter =
710 vtkSmartPointer<vtkTriangleFilter>::New();
711 hullTriangleFilter->SetInputData(hullVTK);
712 hullTriangleFilter->Update();
714 hullVTK = hullTriangleFilter->GetOutput();
716 auto writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
717 writer->SetFileName((fileName +
"_hull.vtp").c_str());
718 writer->SetInputData(hullVTK);
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
float gridDelta
Definition AirGapDeposition.py:21
int counter
Definition Deposition.py:71
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68