TTK
Loading...
Searching...
No Matches
ttkAlgorithm.cpp
Go to the documentation of this file.
1#include <ttkAlgorithm.h>
2#include <ttkMacros.h>
3#include <ttkUtils.h>
4
6#include <Triangulation.h>
8
9#include <vtkCellTypes.h>
10#include <vtkCommand.h>
11#include <vtkDataSet.h>
12
13#ifdef TTK_ENABLE_MPI
15#include <Identifiers.h>
16#include <vtkCellData.h>
17#include <vtkGhostCellsGenerator.h>
18#endif // TTK_ENABLE_MPI
19
20#include <vtkImageData.h>
21#include <vtkInformation.h>
22#include <vtkInformationIntegerKey.h>
23#include <vtkInformationVector.h>
24#include <vtkMultiBlockDataSet.h>
25#include <vtkPointData.h>
26#include <vtkPolyData.h>
27#include <vtkTable.h>
28#include <vtkUnstructuredGrid.h>
29
30#include <vtkCompositeDataPipeline.h>
31
32// Pass input type information key
33#include <vtkInformationKey.h>
34vtkInformationKeyMacro(ttkAlgorithm, SAME_DATA_TYPE_AS_INPUT_PORT, Integer);
35
36// Constructor / Destructor
40
42
43 this->printMsg("Requesting triangulation for '"
44 + std::string(dataSet->GetClassName()) + "'",
46#ifdef TTK_ENABLE_MPI
47 if((ttk::hasInitializedMPI()) && (ttk::isRunningWithMPI())) {
48 if(!hasMPISupport_) {
49 printErr("MPI is not formally supported for this filter :(");
50 printErr("The results are likely to be incorrect.");
51 }
53 }
54#endif // TTK_ENABLE_MPI
55
57 this->debugLevel_, this->CompactTriangulationCacheSize, dataSet);
58
59#ifdef TTK_ENABLE_MPI
60 if(ttk::hasInitializedMPI()) {
61 std::vector<int> tmp{};
62 std::map<int, int> tmpId{};
63 this->MPIPipelinePreconditioning(dataSet, tmp, tmpId, triangulation);
64 this->MPITriangulationPreconditioning(triangulation, dataSet);
65 }
66#endif // TTK_ENABLE_MPI
67
68 if(triangulation)
69 return triangulation;
70
71 this->printErr("Unable to retrieve/initialize triangulation for '"
72 + std::string(dataSet->GetClassName()) + "'");
73
74 return nullptr;
75}
76
77vtkDataArray *ttkAlgorithm::GetOptionalArray(const bool &enforceArrayIndex,
78 const int &arrayIndex,
79 const std::string &arrayName,
80 vtkDataSet *const inputData,
81 const int &inputPort) {
82
83 vtkDataArray *optionalArray = nullptr;
84
85 if(enforceArrayIndex)
86 optionalArray = this->GetInputArrayToProcess(arrayIndex, inputData);
87
88 if(!optionalArray) {
89 this->SetInputArrayToProcess(arrayIndex, inputPort, 0, 0, arrayName.data());
90 optionalArray = this->GetInputArrayToProcess(arrayIndex, inputData);
91 }
92 return optionalArray;
93}
94
95std::string ttkAlgorithm::GetOrderArrayName(vtkDataArray *const array) {
96 return std::string(array->GetName()) + "_Order";
97}
98
99vtkDataArray *
100 ttkAlgorithm::ComputeOrderArray(vtkDataSet *const inputData,
101 vtkDataArray *scalarArray,
102 const int scalarArrayIdx,
103 const bool getGlobalOrder,
104 vtkDataArray *oldOrderArray,
105 ttk::Triangulation *triangulation) {
106
108 auto nVertices = scalarArray->GetNumberOfTuples();
109 if(oldOrderArray != nullptr && getGlobalOrder) {
110 newOrderArray = ttkSimplexIdTypeArray::SafeDownCast(oldOrderArray);
111 } else {
113 newOrderArray->SetName(this->GetOrderArrayName(scalarArray).data());
114 newOrderArray->SetNumberOfComponents(1);
115 newOrderArray->SetNumberOfTuples(nVertices);
116 }
117
118 std::vector<int> neighbors;
119 std::map<int, int> neighborsToId;
120#ifdef TTK_ENABLE_MPI
121 if(ttk::hasInitializedMPI()) {
122 this->MPIGhostPipelinePreconditioning(inputData);
124 inputData, neighbors, neighborsToId, nullptr);
125 }
126 if(ttk::isRunningWithMPI() && getGlobalOrder) {
127 ttk::ArrayPreconditioning arrayPreconditioning
129 arrayPreconditioning.preconditionTriangulation(triangulation);
130 arrayPreconditioning.setGlobalOrder(getGlobalOrder);
131 ttkTypeMacroAT(scalarArray->GetDataType(), triangulation->getType(),
132 (arrayPreconditioning.processScalarArray<T0, T1>(
133 static_cast<const T1 *>(triangulation->getData()),
134 ttkUtils::GetPointer<ttk::SimplexId>(newOrderArray),
135 ttkUtils::GetPointer<T0>(scalarArray), nVertices)));
136 } else {
137 switch(scalarArray->GetDataType()) {
138 vtkTemplateMacro(ttk::preconditionOrderArray(
139 nVertices, static_cast<VTK_TT *>(ttkUtils::GetVoidPointer(scalarArray)),
140 static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(newOrderArray)),
141 this->threadNumber_));
142 }
143 }
144 if(oldOrderArray == nullptr || !getGlobalOrder) {
145 inputData
146 ->GetAttributesAsFieldData(
147 this->GetInputArrayAssociation(scalarArrayIdx, inputData))
148 ->AddArray(newOrderArray);
149 }
150#else
151 switch(scalarArray->GetDataType()) {
152 vtkTemplateMacro(ttk::preconditionOrderArray(
153 nVertices, static_cast<VTK_TT *>(ttkUtils::GetVoidPointer(scalarArray)),
154 static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(newOrderArray)),
155 this->threadNumber_));
156 }
157 inputData
158 ->GetAttributesAsFieldData(
159 this->GetInputArrayAssociation(scalarArrayIdx, inputData))
160 ->AddArray(newOrderArray);
161 TTK_FORCE_USE(triangulation);
162#endif
163 return newOrderArray;
164}
165
167 vtkDataSet *const inputData,
168 vtkDataArray *scalarArray,
169 const int scalarArrayIdx,
170 const bool getGlobalOrder,
171 vtkDataArray *orderArray,
172 ttk::Triangulation *triangulation,
173 const bool enforceOrderArrayIdx) {
174
175 std::string enforcedArray = "";
176 if(enforceOrderArrayIdx) {
177 enforcedArray = " enforced ";
178 }
179#ifdef TTK_ENABLE_MPI
180 if(getGlobalOrder) {
181 if(triangulation->isOrderArrayGlobal(
182 ttkUtils::GetVoidPointer(scalarArray))) {
183 this->printMsg("Retrieved " + enforcedArray + " order array `"
184 + std::string(orderArray->GetName()) + "`.",
186 return orderArray;
187 } else {
188 ttk::Timer timer;
190 this->printWrn("Order array `" + std::string(orderArray->GetName())
191 + "` is local, but a global order array is "
192 "required. Re-computing.");
193
194 this->printMsg("Initializing order array.", 0, 0, this->threadNumber_,
197
198 orderArray
199 = this->ComputeOrderArray(inputData, scalarArray, scalarArrayIdx,
200 getGlobalOrder, orderArray, triangulation);
201
202 triangulation->setIsOrderArrayGlobal(
203 ttkUtils::GetVoidPointer(scalarArray), true);
204
205 this->printMsg("Initializing order array.", 1, timer.getElapsedTime(),
206 this->threadNumber_);
207
209 this->printWrn("TIP: run `ttkArrayPreconditioning` first with "
210 "GlobalOrder enabled");
211 this->printWrn("for improved performances :)");
213 return orderArray;
214 }
215 } else {
216#else
217 TTK_FORCE_USE(inputData);
218 TTK_FORCE_USE(scalarArray);
219 TTK_FORCE_USE(scalarArrayIdx);
220 TTK_FORCE_USE(getGlobalOrder);
221 TTK_FORCE_USE(triangulation);
222#endif // TTK_ENABLE_MPI
223 this->printMsg("Retrieved " + enforcedArray + " order array `"
224 + std::string(orderArray->GetName()) + "`.",
226 return orderArray;
227#ifdef TTK_ENABLE_MPI
228 }
229#endif
230}
231
232vtkDataArray *ttkAlgorithm::GetOrderArray(vtkDataSet *const inputData,
233 const int scalarArrayIdx,
234 ttk::Triangulation *triangulation,
235 const bool getGlobalOrder,
236 const int orderArrayIdx,
237 const bool enforceOrderArrayIdx) {
238
239 auto isValidOrderArray = [](vtkDataArray *const array) {
240 if(!array)
241 return -4;
242
243 if(array->GetNumberOfComponents() != 1)
244 return -3;
245
247 if(array->GetDataType() != temp->GetDataType())
248 return -2;
249
250 const std::string name(array->GetName());
251 if(name.size() < 6 || (name.rfind("_Order") != (name.size() - 6)))
252 return -1;
253
254 return 1;
255 };
256 auto scalarArray = this->GetInputArrayToProcess(scalarArrayIdx, inputData);
257 if(enforceOrderArrayIdx) {
258 auto orderArray = this->GetInputArrayToProcess(orderArrayIdx, inputData);
259 switch(isValidOrderArray(orderArray)) {
260 case -4: {
261 this->printErr("Unable to retrieve enforced order array at idx "
262 + std::to_string(orderArrayIdx) + ".");
263 return nullptr;
264 }
265 case -3: {
266 this->printErr("Retrieved enforced order array `"
267 + std::string(orderArray->GetName())
268 + "` has more than one component.");
269 return nullptr;
270 }
271 case -2: {
272 this->printErr("Enforced order array `"
273 + std::string(orderArray->GetName())
274 + "` is of incorrect type.");
276 this->printErr(" -> use `ttkArrayEditor` to convert data type to `"
277 + std::string(temp->GetDataTypeAsString()) + "`.");
278 return nullptr;
279 }
280 default: {
282 inputData, scalarArray, scalarArrayIdx, getGlobalOrder, orderArray,
283 triangulation, enforceOrderArrayIdx);
284 }
285 }
286 }
287
288 if(!scalarArray) {
289 this->printErr("Unable to retrieve input scalar array for idx "
290 + std::to_string(scalarArrayIdx) + ".");
291 return nullptr;
292 } else if(isValidOrderArray(scalarArray) == 1) {
293 this->printMsg("Retrieved scalar array `"
294 + std::string(scalarArray->GetName())
295 + "` is already an order array.",
297 return scalarArray;
298 }
299
300 auto orderArray = inputData
301 ->GetAttributesAsFieldData(this->GetInputArrayAssociation(
302 scalarArrayIdx, inputData))
303 ->GetArray(this->GetOrderArrayName(scalarArray).data());
304
305 switch(isValidOrderArray(orderArray)) {
306 case -4: {
307 ttk::Timer timer;
309 this->printWrn("No pre-existing order for array:");
310 this->printWrn(" `" + std::string(scalarArray->GetName()) + "`.");
311
312 this->printMsg("Initializing order array.", 0, 0, this->threadNumber_,
315
316 orderArray
317 = this->ComputeOrderArray(inputData, scalarArray, scalarArrayIdx,
318 getGlobalOrder, orderArray, triangulation);
319
320 std::string optionOn = "";
321#ifdef TTK_ENABLE_MPI
322 if(getGlobalOrder) {
323 optionOn = "with GlobalOrder enabled ";
324 }
325 bool isGlobalOrder = getGlobalOrder || (!ttk::isRunningWithMPI());
326 triangulation->setIsOrderArrayGlobal(
327 ttkUtils::GetVoidPointer(scalarArray), isGlobalOrder);
328#endif // TTK_ENABLE_MPI
329 this->printMsg("Initializing order array.", 1, timer.getElapsedTime(),
330 this->threadNumber_);
331
333 this->printWrn("TIP: run `ttkArrayPreconditioning` first");
334 this->printWrn(optionOn + "for improved performances :)");
336
337 return orderArray;
338 }
339
340 case -3: {
341 this->printErr(
342 "Retrieved order array `" + std::string(orderArray->GetName())
343 + "` for scalar array `" + std::string(scalarArray->GetName())
344 + "` has more than one component.");
345 return nullptr;
346 }
347
348 case -2: {
349 this->printErr(
350 "Retrieved order array `" + std::string(orderArray->GetName())
351 + "` for scalar array `" + std::string(scalarArray->GetName())
352 + "` is of incorrect type.");
354 this->printErr(" -> use `ttkArrayEditor` to convert data type to `"
355 + std::string(temp->GetDataTypeAsString()) + "`.");
356 return nullptr;
357 }
358
359 default: {
361 inputData, scalarArray, scalarArrayIdx, getGlobalOrder, orderArray,
362 triangulation, enforceOrderArrayIdx);
363 }
364 }
365}
366
368 ttkAlgorithm::GetIdentifierArrayPtr(const bool &enforceArrayIndex,
369 const int &arrayIndex,
370 const std::string &arrayName,
371 vtkDataSet *const inputData,
372 std::vector<ttk::SimplexId> &spareStorage,
373 const int inputPort,
374 const bool printErr) {
375
376 // fetch data array
377 const auto array = this->GetOptionalArray(
378 enforceArrayIndex, arrayIndex, arrayName, inputData, inputPort);
379 if(array == nullptr) {
380 if(printErr) {
381 this->printErr("Could not find the requested identifiers array");
382 }
383 return {};
384 }
385 if(array->GetNumberOfComponents() != 1) {
386 if(printErr) {
387 this->printErr("Identifiers field must have only one component!");
388 }
389 return {};
390 }
391
392#ifndef TTK_ENABLE_64BIT_IDS
393 if(array->GetDataType() == VTK_ID_TYPE
394 || array->GetDataType() == VTK_LONG_LONG) {
395 this->printMsg(
396 "Converting identifiers field from vtkIdType to SimplexId...");
397 const auto nItems = array->GetNumberOfTuples();
398
399 // fills the vector with the content of the data array converted to
400 // ttk::SimplexId
401 spareStorage.resize(nItems);
402 for(vtkIdType i = 0; i < nItems; ++i) {
403 spareStorage[i] = static_cast<ttk::SimplexId>(array->GetTuple1(i));
404 }
405
406 // return a pointer to the vector internal buffer
407 return spareStorage.data();
408 }
409#else
410 TTK_FORCE_USE(spareStorage);
411#endif // TTK_ENABLE_64BIT_IDS
412
413 // return a pointer to the data array internal buffer
414 return static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(array));
415}
416
417template <class vtkDataType>
418int prepOutput(vtkInformation *info, const std::string &className) {
419 auto output = vtkDataObject::GetData(info);
420 if(!output || !output->IsA(className.data())) {
421 auto newOutput = vtkSmartPointer<vtkDataType>::New();
422 info->Set(vtkDataObject::DATA_OBJECT(), newOutput);
423 }
424 return 1;
425}
426
428 return this->GetOutput(0);
429}
430
431vtkDataSet *ttkAlgorithm::GetOutput(int port) {
432 return vtkDataSet::SafeDownCast(this->GetOutputDataObject(port));
433}
434
435void ttkAlgorithm::SetInputData(vtkDataSet *input) {
436 this->SetInputData(0, input);
437}
438
439void ttkAlgorithm::SetInputData(int index, vtkDataSet *input) {
440 this->SetInputDataInternal(index, input);
441}
442
443void ttkAlgorithm::AddInputData(vtkDataSet *input) {
444 this->AddInputData(0, input);
445}
446
447void ttkAlgorithm::AddInputData(int index, vtkDataSet *input) {
448 this->AddInputDataInternal(index, input);
449}
450
451int ttkAlgorithm::RequestDataObject(vtkInformation *ttkNotUsed(request),
452 vtkInformationVector **inputVector,
453 vtkInformationVector *outputVector) {
454 // for each output
455 for(int i = 0; i < this->GetNumberOfOutputPorts(); ++i) {
456 auto outInfo = outputVector->GetInformationObject(i);
457 if(!outInfo) {
458 this->printErr("Unable to retrieve output vtkDataObject at port "
459 + std::to_string(i));
460 return 0;
461 }
462
463 auto outputPortInfo = this->GetOutputPortInformation(i);
464
465 // always request output type again for dynamic filter outputs
466 if(!this->FillOutputPortInformation(i, outputPortInfo)) {
467 this->printErr("Unable to fill output port information at port "
468 + std::to_string(i));
469 return 0;
470 }
471
472 if(outputPortInfo->Has(ttkAlgorithm::SAME_DATA_TYPE_AS_INPUT_PORT())) {
473 // Set output data type to input data type at specified port
474 auto inPortIndex
475 = outputPortInfo->Get(ttkAlgorithm::SAME_DATA_TYPE_AS_INPUT_PORT());
476 if(inPortIndex < 0 || inPortIndex >= this->GetNumberOfInputPorts()) {
477 this->printErr("Input port index " + std::to_string(inPortIndex)
478 + " specified by 'SAME_DATA_TYPE_AS_INPUT_PORT' key of "
479 "output port is out of range ("
480 + std::to_string(this->GetNumberOfInputPorts())
481 + " input ports).");
482 return 0;
483 }
484 auto inInfo = inputVector[inPortIndex]->GetInformationObject(0);
485 if(!inInfo) {
486 this->printErr(
487 "No information object at port " + std::to_string(inPortIndex)
488 + " specified by 'SAME_DATA_TYPE_AS_INPUT_PORT' key of output port.");
489 return 0;
490 }
491
492 auto input = vtkDataObject::GetData(inInfo);
493 auto output = vtkDataObject::GetData(outInfo);
494
495 if(!output || !output->IsA(input->GetClassName())) {
496 auto newOutput
497 = vtkSmartPointer<vtkDataObject>::Take(input->NewInstance());
498 outputPortInfo->Set(
499 vtkDataObject::DATA_TYPE_NAME(), input->GetClassName());
500 outInfo->Set(vtkDataObject::DATA_OBJECT(), newOutput);
501 }
502 } else {
503 // Explicitly create output by data type name
504 if(!outputPortInfo->Has(vtkDataObject::DATA_TYPE_NAME())) {
505 this->printErr("DATA_TYPE_NAME of output port " + std::to_string(i)
506 + " not specified");
507 return 0;
508 }
509 std::string const outputType
510 = outputPortInfo->Get(vtkDataObject::DATA_TYPE_NAME());
511
512 if(outputType == "vtkUnstructuredGrid") {
513 prepOutput<vtkUnstructuredGrid>(outInfo, outputType);
514 } else if(outputType == "vtkPolyData") {
515 prepOutput<vtkPolyData>(outInfo, outputType);
516 } else if(outputType == "vtkMultiBlockDataSet") {
517 prepOutput<vtkMultiBlockDataSet>(outInfo, outputType);
518 } else if(outputType == "vtkTable") {
519 prepOutput<vtkTable>(outInfo, outputType);
520 } else if(outputType == "vtkImageData") {
521 prepOutput<vtkImageData>(outInfo, outputType);
522 } else {
523 this->printErr("Unsupported data type for output[" + std::to_string(i)
524 + "]: " + outputType);
525 return 0;
526 }
527 }
528
529 this->printMsg(
530 "Created '"
531 + std::string(outputPortInfo->Get(vtkDataObject::DATA_TYPE_NAME()))
532 + "' at output port " + std::to_string(i),
534 }
535
536 return 1;
537}
538
539#ifdef TTK_ENABLE_MPI
540
541int ttkAlgorithm::updateMPICommunicator(vtkDataSet *input) {
542 if(input == nullptr) {
543 return 0;
544 }
545 int isEmpty
546 = input->GetNumberOfCells() == 0 || input->GetNumberOfPoints() == 0;
547 int oldSize = ttk::MPIsize_;
548 int oldRank = ttk::MPIrank_;
549 MPI_Comm_split(MPI_COMM_WORLD, isEmpty, 0, &ttk::MPIcomm_);
550 MPI_Comm_rank(ttk::MPIcomm_, &ttk::MPIrank_);
551 MPI_Comm_size(ttk::MPIcomm_, &ttk::MPIsize_);
552 if(oldSize != ttk::MPIsize_) {
553 std::vector<int> newToOldRanks(ttk::MPIsize_);
554 MPI_Allgather(&oldRank, 1, MPI_INTEGER, newToOldRanks.data(), 1,
555 MPI_INTEGER, ttk::MPIcomm_);
556 std::map<int, int> oldToNewRanks;
557 for(int i = 0; i < ttk::MPIsize_; i++) {
558 oldToNewRanks[newToOldRanks[i]] = i;
559 }
560 int *vertexRankArray
561 = ttkUtils::GetPointer<int>(input->GetPointData()->GetArray("RankArray"));
562 if(vertexRankArray != nullptr) {
563 for(int i = 0; i < input->GetNumberOfPoints(); i++) {
564 vertexRankArray[i] = oldToNewRanks[vertexRankArray[i]];
565 }
566 }
567 int *cellRankArray
568 = ttkUtils::GetPointer<int>(input->GetCellData()->GetArray("RankArray"));
569 if(cellRankArray != nullptr) {
570 for(int i = 0; i < input->GetNumberOfCells(); i++) {
571 cellRankArray[i] = oldToNewRanks[cellRankArray[i]];
572 }
573 }
574 }
576 return isEmpty;
577}
578
580 ttk::SimplexId simplexNumber,
581 unsigned char *ghost,
582 int *rankArray) {
583 ttk::SimplexId ghostNumber = 0;
584 if(rankArray != nullptr) {
585#ifdef TTK_ENABLE_OPENMP
586#pragma omp parallel for reduction(+ : ghostNumber)
587#endif // TTK_ENABLE_OPENMP
588 for(ttk::SimplexId i = 0; i < simplexNumber; i++) {
589 if(rankArray[i] != ttk::MPIrank_) {
590 ghostNumber++;
591 }
592 }
593 } else {
594 if(ghost != nullptr) {
595#ifdef TTK_ENABLE_OPENMP
596#pragma omp parallel for reduction(+ : ghostNumber)
597#endif // TTK_ENABLE_OPENMP
598 for(ttk::SimplexId i = 0; i < simplexNumber; i++) {
599 if(ghost[i] == 1) {
600 ghostNumber++;
601 }
602 }
603 }
604 }
605
606 ttk::SimplexId realSimplexNumber = simplexNumber - ghostNumber;
607 auto minmax = std::minmax_element(globalIds, globalIds + simplexNumber);
608 ttk::LongSimplexId min = globalIds[minmax.first - globalIds];
609 ttk::LongSimplexId max = globalIds[minmax.second - globalIds];
610 ttk::SimplexId globalSimplexNumber;
611 ttk::LongSimplexId globalMin;
612 ttk::LongSimplexId globalMax;
613 MPI_Allreduce(&realSimplexNumber, &globalSimplexNumber, 1,
614 ttk::getMPIType(realSimplexNumber), MPI_SUM, ttk::MPIcomm_);
615 MPI_Allreduce(
616 &min, &globalMin, 1, ttk::getMPIType(min), MPI_MIN, ttk::MPIcomm_);
617 MPI_Allreduce(
618 &max, &globalMax, 1, ttk::getMPIType(max), MPI_MAX, ttk::MPIcomm_);
619
620 return (globalSimplexNumber == globalMax + 1 && globalMin == 0);
621};
622
624 vtkDataSet *input,
625 std::unordered_map<ttk::SimplexId, ttk::SimplexId> &vertGtoL,
626 std::vector<int> &neighborRanks,
627 std::map<int, int> &neighborsToId) {
628
629 ttk::Identifiers identifiers;
630
631 vtkNew<vtkIdTypeArray> vtkVertexIdentifiers{};
632
633 vtkNew<vtkIdTypeArray> vtkCellIdentifiers{};
634 vtkVertexIdentifiers->SetName("GlobalPointIds");
635 vtkVertexIdentifiers->SetNumberOfComponents(1);
636 vtkVertexIdentifiers->SetNumberOfTuples(input->GetNumberOfPoints());
637 vtkVertexIdentifiers->Fill(-1);
638
639 identifiers.setVertexIdentifiers(
640 ttkUtils::GetPointer<ttk::LongSimplexId>(vtkVertexIdentifiers));
641
642 vtkCellIdentifiers->SetName("GlobalCellIds");
643 vtkCellIdentifiers->SetNumberOfComponents(1);
644 vtkCellIdentifiers->SetNumberOfTuples(input->GetNumberOfCells());
645 vtkCellIdentifiers->Fill(-1);
646
647 identifiers.setCellIdentifiers(
648 ttkUtils::GetPointer<ttk::LongSimplexId>(vtkCellIdentifiers));
649
650 int vertexNumber = input->GetNumberOfPoints();
651 identifiers.setVertexNumber(vertexNumber);
652 int cellNumber = input->GetNumberOfCells();
653 identifiers.setCellNumber(cellNumber);
654 int status = 0;
655
656 double *boundingBox = input->GetBounds();
657 identifiers.setBounds(boundingBox);
658 identifiers.initializeNeighbors(boundingBox, neighborRanks, neighborsToId);
659 if(ttk::isRunningWithMPI()) {
660 switch(input->GetDataObjectType()) {
661 case VTK_UNSTRUCTURED_GRID:
662 case VTK_POLY_DATA: {
663
664 identifiers.setOutdatedGlobalPointIds(
665 ttkUtils::GetPointer<ttk::LongSimplexId>(
666 input->GetPointData()->GetGlobalIds()));
667 identifiers.setOutdatedGlobalCellIds(
668 ttkUtils::GetPointer<ttk::LongSimplexId>(
669 input->GetCellData()->GetGlobalIds()));
670 identifiers.setVertexRankArray(ttkUtils::GetPointer<int>(
671 input->GetPointData()->GetArray("RankArray")));
672 int *cellRankArray = ttkUtils::GetPointer<int>(
673 input->GetCellData()->GetArray("RankArray"));
674 identifiers.setCellRankArray(cellRankArray);
675 identifiers.setVertGhost(ttkUtils::GetPointer<unsigned char>(
676 input->GetPointData()->GetArray("vtkGhostType")));
677 unsigned char *cellGhost = ttkUtils::GetPointer<unsigned char>(
678 input->GetCellData()->GetArray("vtkGhostType"));
679 identifiers.setCellGhost(cellGhost);
680 vtkPointSet *pointSet = vtkPointSet::SafeDownCast(input);
681 identifiers.setPointSet(static_cast<float *>(
682 ttkUtils::GetVoidPointer(pointSet->GetPoints())));
683 vtkCellArray *cells = nullptr;
684 switch(input->GetDataObjectType()) {
685 case VTK_UNSTRUCTURED_GRID: {
686 auto dataSetAsUG = vtkUnstructuredGrid::SafeDownCast(input);
687 cells = dataSetAsUG->GetCells();
688 break;
689 }
690 case VTK_POLY_DATA: {
691 auto dataSetAsPD = vtkPolyData::SafeDownCast(input);
692 cells
693 = dataSetAsPD->GetNumberOfPolys() > 0 ? dataSetAsPD->GetPolys()
694 : dataSetAsPD->GetNumberOfLines() > 0 ? dataSetAsPD->GetLines()
695 : dataSetAsPD->GetVerts();
696 break;
697 }
698 default: {
699 this->printErr("Unable to get cells for `"
700 + std::string(input->GetClassName()) + "`");
701 }
702 }
703 if(cells == nullptr) {
704 return 0;
705 }
706 if(!cells->IsStorage64Bit()) {
707 if(cells->CanConvertTo64BitStorage()) {
708 this->printWrn("Converting the cell array to 64-bit storage");
709 bool success = cells->ConvertTo64BitStorage();
710 if(!success) {
711 this->printErr(
712 "Error converting the provided cell array to 64-bit storage");
713 return -1;
714 }
715 } else {
716 this->printErr(
717 "Cannot convert the provided cell array to 64-bit storage");
718 return -1;
719 }
720 }
721
722 identifiers.setConnectivity(ttkUtils::GetPointer<ttk::LongSimplexId>(
723 cells->GetConnectivityArray()));
724
725 std::vector<std::vector<ttk::SimplexId>> pointsToCells(vertexNumber);
726 vtkIdList *cellList = vtkIdList::New();
727 if(cellRankArray != nullptr) {
728 for(ttk::SimplexId i = 0; i < vertexNumber; i++) {
729 input->GetPointCells(i, cellList);
730 for(int j = 0; j < cellList->GetNumberOfIds(); j++) {
731 if(cellRankArray[cellList->GetId(j)] == ttk::MPIrank_) {
732 pointsToCells[i].push_back(cellList->GetId(j));
733 }
734 }
735 }
736 } else {
737 for(ttk::SimplexId i = 0; i < vertexNumber; i++) {
738 input->GetPointCells(i, cellList);
739 for(int j = 0; j < cellList->GetNumberOfIds(); j++) {
740 if(cellGhost[cellList->GetId(j)] == 0) {
741 pointsToCells[i].push_back(cellList->GetId(j));
742 }
743 }
744 }
745 }
746 identifiers.setPointsToCells(pointsToCells);
747
748 identifiers.initializeMPITypes();
749 identifiers.setVertGtoL(&vertGtoL);
750 vtkIdList *pointCell = vtkIdList::New();
751 input->GetCellPoints(0, pointCell);
752 int nbPoints = pointCell->GetNumberOfIds();
753 identifiers.setDomainDimension(nbPoints - 1);
754 identifiers.buildKDTree();
755 status = identifiers.executePolyData();
756 break;
757 }
758 case VTK_IMAGE_DATA: {
759 vtkImageData *data = vtkImageData::SafeDownCast(input);
760 identifiers.setDims(data->GetDimensions());
761 identifiers.setSpacing(data->GetSpacing());
762 status = identifiers.executeImageData();
763 break;
764 }
765 default: {
766 this->printErr("Unable to triangulate `"
767 + std::string(input->GetClassName()) + "`");
768 }
769 }
770 } else {
771 status = identifiers.executeSequential();
772 }
773
774 if(status < 1) {
775 printErr("Global identifier generation failed");
776 return -1;
777 }
778
779 // Add VTK objects to the data set
780
781 input->GetPointData()->SetGlobalIds(vtkVertexIdentifiers);
782 input->GetCellData()->SetGlobalIds(vtkCellIdentifiers);
783 return 0;
784}
785
786void ttkAlgorithm::MPIGhostPipelinePreconditioning(vtkDataSet *input) {
787
788 vtkNew<vtkGhostCellsGenerator> generator;
789 if(ttk::isRunningWithMPI()
790 && (!input->HasAnyGhostCells()
791 && ((input->GetPointData()->GetArray("RankArray") == nullptr)
792 || (input->GetCellData()->GetArray("RankArray") == nullptr)))) {
793 generator->SetInputData(input);
794 generator->BuildIfRequiredOff();
795 generator->SetNumberOfGhostLayers(1);
796 generator->Update();
797 input->ShallowCopy(generator->GetOutputDataObject(0));
798 input->GetPointData()->AddArray(
799 generator->GetOutputDataObject(0)->GetGhostArray(0));
800 input->GetCellData()->AddArray(
801 generator->GetOutputDataObject(0)->GetGhostArray(1));
802 }
803}
804
806 vtkDataSet *input,
807 std::vector<int> &neighbors,
808 std::map<int, int> &neighToId,
809 ttk::Triangulation *triangulation) {
810
811 ttk::SimplexId vertexNumber = input->GetNumberOfPoints();
812 ttk::SimplexId cellNumber = input->GetNumberOfCells();
813
814 if((input->GetDataObjectType() == VTK_POLY_DATA
815 || input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)) {
816
817 if((ttk::hasInitializedMPI()) && (ttk::isRunningWithMPI())) {
819 printWrn("The distribution by VTK of Unstructured");
820 printWrn("Grids and Poly Data has been reported");
821 printWrn("to be affected by bugs (at least up");
822 printWrn("to ParaView 5.10.1).");
823
824 if((input->GetCellData()->GetGlobalIds() == nullptr)
825 || (input->GetPointData()->GetGlobalIds() == nullptr)) {
826
827 printWrn("=> Global identifiers may be incorrect.");
828 }
829 if((input->GetPointData()->GetArray("RankArray") == nullptr)
830 || (input->GetCellData()->GetArray("RankArray") == nullptr)) {
831 printWrn("=> Rank arrays may be incorrect.");
832 }
834 }
835 }
836
837 // Get the neighbor ranks
838 std::vector<int> &neighborRanks{
839 triangulation != nullptr ? triangulation->getNeighborRanks() : neighbors};
840 std::map<int, int> &neighborsToId{
841 triangulation != nullptr ? triangulation->getNeighborsToId() : neighToId};
842
843 double *boundingBox = input->GetBounds();
844 if(triangulation != nullptr) {
845 triangulation->createMetaGrid(boundingBox);
846 }
847
848 if(neighborRanks.empty()) {
849 ttk::preconditionNeighborsUsingBoundingBox(
850 boundingBox, neighborRanks, neighborsToId);
851 }
852
853 // Checks if global ids are valid
854 ttk::LongSimplexId *globalPointIds = ttkUtils::GetPointer<ttk::LongSimplexId>(
855 input->GetPointData()->GetGlobalIds());
856 ttk::LongSimplexId *globalCellIds = ttkUtils::GetPointer<ttk::LongSimplexId>(
857 input->GetCellData()->GetGlobalIds());
858
859 bool pointValidity{false};
860 bool cellValidity{false};
861 if((triangulation != nullptr
862 && (triangulation->getType() == ttk::Triangulation::Type::EXPLICIT
863 || triangulation->getType() == ttk::Triangulation::Type::COMPACT))
864 || triangulation == nullptr) {
865 if(globalPointIds != nullptr) {
866 unsigned char *ghostPoints = ttkUtils::GetPointer<unsigned char>(
867 input->GetPointData()->GetArray("vtkGhostType"));
868 int *vertexRankArray = ttkUtils::GetPointer<int>(
869 input->GetPointData()->GetArray("RankArray"));
870 pointValidity = checkGlobalIdValidity(
871 globalPointIds, vertexNumber, ghostPoints, vertexRankArray);
872 }
873 if(pointValidity && globalCellIds != nullptr) {
874
875 unsigned char *ghostCells = ttkUtils::GetPointer<unsigned char>(
876 input->GetCellData()->GetArray("vtkGhostType"));
877 int *cellRankArray = ttkUtils::GetPointer<int>(
878 input->GetCellData()->GetArray("RankArray"));
879 cellValidity = checkGlobalIdValidity(
880 globalCellIds, cellNumber, ghostCells, cellRankArray);
881 }
882 } else {
883 pointValidity = true;
884 cellValidity = true;
885 }
886
887 // If the global ids are not valid, they are computed again
888 if(!pointValidity || !cellValidity) {
889 if(triangulation != nullptr) {
890 if(triangulation->getType() == ttk::Triangulation::Type::EXPLICIT
891 || triangulation->getType() == ttk::Triangulation::Type::COMPACT) {
892 this->GenerateGlobalIds(input, triangulation->getVertexGlobalIdMap(),
893 neighborRanks, neighborsToId);
894 }
895 } else {
896 std::unordered_map<ttk::SimplexId, ttk::SimplexId> vertGtoL{};
897 this->GenerateGlobalIds(input, vertGtoL, neighborRanks, neighborsToId);
898 }
899 }
900}
901
903 ttk::Triangulation *triangulation, vtkDataSet *input) {
904
905 const auto pd{input->GetPointData()};
906 if(pd == nullptr) {
907 triangulation->printWrn("No point data on input object");
908 } else {
909 // provide "GlobalPointIds" & "vtkGhostType" point data arrays
910 // to the triangulation
911 triangulation->setVertsGlobalIds(
912 ttkUtils::GetPointer<ttk::LongSimplexId>(pd->GetGlobalIds()));
913 triangulation->setVertexGhostArray(
914 ttkUtils::GetPointer<unsigned char>(pd->GetArray("vtkGhostType")));
915 int *vertexRankArray = ttkUtils::GetPointer<int>(pd->GetArray("RankArray"));
916 if(vertexRankArray != nullptr) {
917 triangulation->setVertexRankArray(vertexRankArray);
918 }
919 triangulation->preconditionDistributedVertices();
920 }
921
922 const auto cd{input->GetCellData()};
923 if(cd == nullptr) {
924 triangulation->printWrn("No cell data on input object");
925 } else {
926 // provide "GlobalCellIds" & "vtkGhostType" cell data arrays to
927 // the triangulation
928 triangulation->setCellsGlobalIds(
929 ttkUtils::GetPointer<ttk::LongSimplexId>(cd->GetGlobalIds()));
930 triangulation->setCellGhostArray(
931 ttkUtils::GetPointer<unsigned char>(cd->GetArray("vtkGhostType")));
932 int *cellRankArray = ttkUtils::GetPointer<int>(cd->GetArray("RankArray"));
933 if(cellRankArray != nullptr) {
934 triangulation->setCellRankArray(cellRankArray);
935 }
936 triangulation->preconditionDistributedCells();
937 }
938}
939
940#endif // TTK_ENABLE_MPI
941
942//==============================================================================
943int ttkAlgorithm::ProcessRequest(vtkInformation *request,
944 vtkInformationVector **inputVector,
945 vtkInformationVector *outputVector) {
946 // 1. Pass
947 if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA_OBJECT())) {
948 this->printMsg(
949 "Processing REQUEST_DATA_OBJECT", ttk::debug::Priority::VERBOSE);
950 return this->RequestDataObject(request, inputVector, outputVector);
951 }
952
953 // 2. Pass
954 if(request->Has(vtkCompositeDataPipeline::REQUEST_INFORMATION())) {
955 this->printMsg(
956 "Processing REQUEST_INFORMATION", ttk::debug::Priority::VERBOSE);
957 return this->RequestInformation(request, inputVector, outputVector);
958 }
959
960 // 3. Pass
961 if(request->Has(vtkCompositeDataPipeline::REQUEST_UPDATE_TIME())) {
962 this->printMsg(
963 "Processing REQUEST_UPDATE_TIME", ttk::debug::Priority::VERBOSE);
964 return this->RequestUpdateTime(request, inputVector, outputVector);
965 }
966
967 // 4. Pass
968 if(request->Has(
969 vtkCompositeDataPipeline::REQUEST_TIME_DEPENDENT_INFORMATION())) {
970 this->printMsg("Processing REQUEST_TIME_DEPENDENT_INFORMATION",
973 request, inputVector, outputVector);
974 }
975
976 // 5. Pass
977 if(request->Has(vtkCompositeDataPipeline::REQUEST_UPDATE_EXTENT())) {
978 this->printMsg(
979 "Processing REQUEST_UPDATE_EXTENT", ttk::debug::Priority::VERBOSE);
980 return this->RequestUpdateExtent(request, inputVector, outputVector);
981 }
982
983 // 6. Pass
984 if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA_NOT_GENERATED())) {
985 this->printMsg(
986 "Processing REQUEST_DATA_NOT_GENERATED", ttk::debug::Priority::VERBOSE);
987 return this->RequestDataNotGenerated(request, inputVector, outputVector);
988 }
989
990 // 7. Pass
991 if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA())) {
992 this->printMsg("Processing REQUEST_DATA", ttk::debug::Priority::VERBOSE);
994#ifdef TTK_ENABLE_MPI
995 if(ttk::hasInitializedMPI() && inputVector != nullptr) {
996 if(this->updateMPICommunicator(vtkDataSet::GetData(inputVector[0], 0))) {
997 return 1;
998 };
999 }
1000#endif // TTK_ENABLE_MPI
1001 return this->RequestData(request, inputVector, outputVector);
1002 }
1003
1004 this->printErr("Unsupported pipeline pass:");
1005 request->Print(cout);
1006
1007 return 0;
1008}
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition BaseClass.h:57
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition BaseClass.h:47
Baseclass of all VTK filters that wrap ttk modules.
void SetInputData(vtkDataSet *)
static vtkInformationIntegerKey * SAME_DATA_TYPE_AS_INPUT_PORT()
bool checkGlobalIdValidity(ttk::LongSimplexId *globalIds, ttk::SimplexId simplexNumber, unsigned char *ghost, int *rankArray)
ttk::SimplexId * GetIdentifierArrayPtr(const bool &enforceArrayIndex, const int &arrayIndex, const std::string &arrayName, vtkDataSet *const inputData, std::vector< ttk::SimplexId > &spareStorage, const int inputPort=0, const bool printErr=true)
void AddInputData(vtkDataSet *)
virtual int RequestUpdateTime(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
~ttkAlgorithm() override
virtual int RequestDataNotGenerated(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
float CompactTriangulationCacheSize
virtual int RequestData(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
void MPIGhostPipelinePreconditioning(vtkDataSet *input)
virtual int RequestInformation(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
int ProcessRequest(vtkInformation *request, vtkInformationVector **inputVectors, vtkInformationVector *outputVector) override
int GenerateGlobalIds(vtkDataSet *input, std::unordered_map< ttk::SimplexId, ttk::SimplexId > &vertGtoL, std::vector< int > &neighborRanks, std::map< int, int > &neighborsToId)
vtkDataArray * checkForGlobalAndComputeOrderArray(vtkDataSet *const inputData, vtkDataArray *scalarArray, const int scalarArrayIdx, const bool getGlobalOrder, vtkDataArray *orderArray, ttk::Triangulation *triangulation, const bool enforceOrderArrayIdx)
virtual int RequestDataObject(vtkInformation *request, vtkInformationVector **inputVectors, vtkInformationVector *outputVector)
virtual int RequestUpdateTimeDependentInformation(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
ttk::Triangulation * GetTriangulation(vtkDataSet *dataSet)
virtual int RequestUpdateExtent(vtkInformation *ttkNotUsed(request), vtkInformationVector **ttkNotUsed(inputVectors), vtkInformationVector *ttkNotUsed(outputVector))
void MPITriangulationPreconditioning(ttk::Triangulation *triangulation, vtkDataSet *input)
vtkDataSet * GetOutput()
vtkDataArray * GetOrderArray(vtkDataSet *const inputData, const int scalarArrayIdx, ttk::Triangulation *triangulation, const bool getGlobalOrder=false, const int orderArrayIdx=0, const bool enforceOrderArrayIdx=false)
void MPIPipelinePreconditioning(vtkDataSet *input, std::vector< int > &neighbors, std::map< int, int > &neighToId, ttk::Triangulation *triangulation=nullptr)
int FillOutputPortInformation(int ttkNotUsed(port), vtkInformation *ttkNotUsed(info)) override
vtkDataArray * ComputeOrderArray(vtkDataSet *const inputData, vtkDataArray *scalarArray, const int scalarArrayIdx, const bool getGlobalOrder, vtkDataArray *oldOrderArray, ttk::Triangulation *triangulation)
vtkDataArray * GetOptionalArray(const bool &enforceArrayIndex, const int &arrayIndex, const std::string &arrayName, vtkDataSet *const inputData, const int &inputPort=0)
static std::string GetOrderArrayName(vtkDataArray *const array)
static ttk::Triangulation * GetTriangulation(int debugLevel, float cacheRatio, vtkDataSet *object)
static void * GetVoidPointer(vtkDataArray *array, vtkIdType start=0)
Definition ttkUtils.cpp:226
int preconditionTriangulation(AbstractTriangulation *triangulation)
int processScalarArray(const triangulationType *triangulation, ttk::SimplexId *orderArray, const DT *scalarArray, const size_t nVerts) const
int debugLevel_
Definition Debug.h:379
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:159
std::string debugMsgNamePrefix_
Definition Debug.h:384
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
void setCellIdentifiers(ttk::LongSimplexId *cellIdentifiers)
Definition Identifiers.h:97
void setVertexNumber(const SimplexId &vertexNumber)
Definition Identifiers.h:85
void setCellNumber(const SimplexId &cellNumber)
Definition Identifiers.h:89
int executeSequential()
Generates global ids for all data set type in sequential.
void setVertexIdentifiers(ttk::LongSimplexId *vertexIdentifiers)
Definition Identifiers.h:93
double getElapsedTime()
Definition Timer.h:15
Triangulation is a class that provides time and memory efficient traversal methods on triangulations ...
AbstractTriangulation * getData()
Triangulation::Type getType() const
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
void preconditionOrderArray(const size_t nVerts, const scalarType *const scalars, SimplexId *const order, const int nThreads=ttk::globalThreadNumber_)
Precondition an order array to be consumed by the base layer API.
long long int LongSimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:15
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
vtkStandardNewMacro(ttkAlgorithm)
vtkInformationKeyMacro(ttkAlgorithm, SAME_DATA_TYPE_AS_INPUT_PORT, Integer)
int prepOutput(vtkInformation *info, const std::string &className)
#define ttkTypeMacroAT(group0, group1, call)
Definition ttkMacros.h:272
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)