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