ttkAlgorithm.cpp

Go to the documentation of this file.

#include <ttkAlgorithm.h>
#include <ttkMacros.h>
#include <ttkUtils.h>
 
#include <OrderDisambiguation.h>
#include <Triangulation.h>
#include <ttkTriangulationFactory.h>
 
#include <vtkCellTypes.h>
#include <vtkCommand.h>
#include <vtkDataSet.h>
 
#ifdef TTK_ENABLE_MPI
#include <ArrayPreconditioning.h>
#include <Identifiers.h>
#include <vtkCellData.h>
#include <vtkGhostCellsGenerator.h>
#endif // TTK_ENABLE_MPI
 
#include <vtkImageData.h>
#include <vtkInformation.h>
#include <vtkInformationIntegerKey.h>
#include <vtkInformationVector.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkTable.h>
#include <vtkUnstructuredGrid.h>
 
#include <vtkCompositeDataPipeline.h>
 
using namespace std;
 
// Pass input type information key
#include <vtkInformationKey.h>
vtkInformationKeyMacro(ttkAlgorithm, SAME_DATA_TYPE_AS_INPUT_PORT, Integer);
 
// Constructor / Destructor
vtkStandardNewMacro(ttkAlgorithm);
ttkAlgorithm::ttkAlgorithm() = default;
ttkAlgorithm::~ttkAlgorithm() = default;
 
ttk::Triangulation *ttkAlgorithm::GetTriangulation(vtkDataSet *dataSet) {
 
  this->printMsg("Requesting triangulation for '"
                   + std::string(dataSet->GetClassName()) + "'",
                 ttk::debug::Priority::DETAIL);
#ifdef TTK_ENABLE_MPI
  if((ttk::hasInitializedMPI()) && (ttk::isRunningWithMPI())) {
    if(!hasMPISupport_) {
      printErr("MPI is not formally supported for this filter :(");
      printErr("The results are likely to be incorrect.");
    }
    this->MPIGhostPipelinePreconditioning(dataSet);
  }
#endif // TTK_ENABLE_MPI
 
  auto triangulation = ttkTriangulationFactory::GetTriangulation(
    this->debugLevel_, this->CompactTriangulationCacheSize, dataSet);
 
#ifdef TTK_ENABLE_MPI
  if(ttk::hasInitializedMPI()) {
    std::vector<int> tmp{};
    std::map<int, int> tmpId{};
    this->MPIPipelinePreconditioning(dataSet, tmp, tmpId, triangulation);
    this->MPITriangulationPreconditioning(triangulation, dataSet);
  }
#endif // TTK_ENABLE_MPI
 
  if(triangulation)
    return triangulation;
 
  this->printErr("Unable to retrieve/initialize triangulation for '"
                 + std::string(dataSet->GetClassName()) + "'");
 
  return nullptr;
}
 
vtkDataArray *ttkAlgorithm::GetOptionalArray(const bool &enforceArrayIndex,
                                             const int &arrayIndex,
                                             const std::string &arrayName,
                                             vtkDataSet *const inputData,
                                             const int &inputPort) {
 
  vtkDataArray *optionalArray = nullptr;
 
  if(enforceArrayIndex)
    optionalArray = this->GetInputArrayToProcess(arrayIndex, inputData);
 
  if(!optionalArray) {
    this->SetInputArrayToProcess(arrayIndex, inputPort, 0, 0, arrayName.data());
    optionalArray = this->GetInputArrayToProcess(arrayIndex, inputData);
  }
  return optionalArray;
}
 
std::string ttkAlgorithm::GetOrderArrayName(vtkDataArray *const array) {
  return std::string(array->GetName()) + "_Order";
}
 
vtkDataArray *
  ttkAlgorithm::ComputeOrderArray(vtkDataSet *const inputData,
                                  vtkDataArray *scalarArray,
                                  const int scalarArrayIdx,
                                  const bool getGlobalOrder,
                                  vtkDataArray *oldOrderArray,
                                  ttk::Triangulation *triangulation) {
 
  vtkSmartPointer<ttkSimplexIdTypeArray> newOrderArray;
  auto nVertices = scalarArray->GetNumberOfTuples();
  if(oldOrderArray != nullptr && getGlobalOrder) {
    newOrderArray = ttkSimplexIdTypeArray::SafeDownCast(oldOrderArray);
  } else {
    newOrderArray = vtkSmartPointer<ttkSimplexIdTypeArray>::New();
    newOrderArray->SetName(this->GetOrderArrayName(scalarArray).data());
    newOrderArray->SetNumberOfComponents(1);
    newOrderArray->SetNumberOfTuples(nVertices);
  }
 
  std::vector<int> neighbors;
  std::map<int, int> neighborsToId;
#ifdef TTK_ENABLE_MPI
  if(ttk::hasInitializedMPI()) {
    this->MPIGhostPipelinePreconditioning(inputData);
    this->MPIPipelinePreconditioning(
      inputData, neighbors, neighborsToId, nullptr);
  }
  if(ttk::isRunningWithMPI() && getGlobalOrder) {
    ttk::ArrayPreconditioning arrayPreconditioning
      = ttk::ArrayPreconditioning();
    arrayPreconditioning.preconditionTriangulation(triangulation);
    arrayPreconditioning.setGlobalOrder(getGlobalOrder);
    ttkTypeMacroAT(scalarArray->GetDataType(), triangulation->getType(),
                   (arrayPreconditioning.processScalarArray<T0, T1>(
                     static_cast<const T1 *>(triangulation->getData()),
                     ttkUtils::GetPointer<ttk::SimplexId>(newOrderArray),
                     ttkUtils::GetPointer<T0>(scalarArray), nVertices)));
  } else {
    switch(scalarArray->GetDataType()) {
      vtkTemplateMacro(ttk::preconditionOrderArray(
        nVertices, static_cast<VTK_TT *>(ttkUtils::GetVoidPointer(scalarArray)),
        static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(newOrderArray)),
        this->threadNumber_));
    }
  }
  if(oldOrderArray == nullptr || !getGlobalOrder) {
    inputData
      ->GetAttributesAsFieldData(
        this->GetInputArrayAssociation(scalarArrayIdx, inputData))
      ->AddArray(newOrderArray);
  }
#else
 
  switch(scalarArray->GetDataType()) {
    vtkTemplateMacro(ttk::preconditionOrderArray(
      nVertices, static_cast<VTK_TT *>(ttkUtils::GetVoidPointer(scalarArray)),
      static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(newOrderArray)),
      this->threadNumber_));
  }
 
  inputData
    ->GetAttributesAsFieldData(
      this->GetInputArrayAssociation(scalarArrayIdx, inputData))
    ->AddArray(newOrderArray);
  TTK_FORCE_USE(triangulation);
#endif
  return newOrderArray;
}
 
vtkDataArray *ttkAlgorithm::checkForGlobalAndComputeOrderArray(
  vtkDataSet *const inputData,
  vtkDataArray *scalarArray,
  const int scalarArrayIdx,
  const bool getGlobalOrder,
  vtkDataArray *orderArray,
  ttk::Triangulation *triangulation,
  const bool enforceOrderArrayIdx) {
 
  std::string enforcedArray = "";
  if(enforceOrderArrayIdx) {
    enforcedArray = " enforced ";
  }
#ifdef TTK_ENABLE_MPI
  if(getGlobalOrder) {
    if(triangulation->isOrderArrayGlobal(
         ttkUtils::GetVoidPointer(scalarArray))) {
      this->printMsg("Retrieved " + enforcedArray + " order array `"
                       + std::string(orderArray->GetName()) + "`.",
                     ttk::debug::Priority::DETAIL);
      return orderArray;
    } else {
      ttk::Timer timer;
      printMsg(ttk::debug::Separator::L2);
      this->printWrn("Order array `" + std::string(orderArray->GetName())
                     + "` is local, but a global order array is "
                       "required. Re-computing.");
 
      this->printMsg("Initializing order array.", 0, 0, this->threadNumber_,
                     ttk::debug::LineMode::REPLACE);
      printMsg(ttk::debug::Separator::L2);
 
      orderArray
        = this->ComputeOrderArray(inputData, scalarArray, scalarArrayIdx,
                                  getGlobalOrder, orderArray, triangulation);
 
      triangulation->setIsOrderArrayGlobal(
        ttkUtils::GetVoidPointer(scalarArray), true);
 
      this->printMsg("Initializing order array.", 1, timer.getElapsedTime(),
                     this->threadNumber_);
 
      printMsg(ttk::debug::Separator::L2);
      this->printWrn("TIP: run `ttkArrayPreconditioning` first with "
                     "GlobalOrder enabled");
      this->printWrn("for improved performances :)");
      printMsg(ttk::debug::Separator::L2);
      return orderArray;
    }
  } else {
#else
  TTK_FORCE_USE(inputData);
  TTK_FORCE_USE(scalarArray);
  TTK_FORCE_USE(scalarArrayIdx);
  TTK_FORCE_USE(getGlobalOrder);
  TTK_FORCE_USE(triangulation);
#endif // TTK_ENABLE_MPI
    this->printMsg("Retrieved " + enforcedArray + " order array `"
                     + std::string(orderArray->GetName()) + "`.",
                   ttk::debug::Priority::DETAIL);
    return orderArray;
#ifdef TTK_ENABLE_MPI
  }
#endif
}
 
vtkDataArray *ttkAlgorithm::GetOrderArray(vtkDataSet *const inputData,
                                          const int scalarArrayIdx,
                                          ttk::Triangulation *triangulation,
                                          const bool getGlobalOrder,
                                          const int orderArrayIdx,
                                          const bool enforceOrderArrayIdx) {
 
  auto isValidOrderArray = [](vtkDataArray *const array) {
    if(!array)
      return -4;
 
    if(array->GetNumberOfComponents() != 1)
      return -3;
 
    auto temp = vtkSmartPointer<ttkSimplexIdTypeArray>::New();
    if(array->GetDataType() != temp->GetDataType())
      return -2;
 
    const std::string name(array->GetName());
    if(name.size() < 6 || (name.rfind("_Order") != (name.size() - 6)))
      return -1;
 
    return 1;
  };
  auto scalarArray = this->GetInputArrayToProcess(scalarArrayIdx, inputData);
  if(enforceOrderArrayIdx) {
    auto orderArray = this->GetInputArrayToProcess(orderArrayIdx, inputData);
    switch(isValidOrderArray(orderArray)) {
      case -4: {
        this->printErr("Unable to retrieve enforced order array at idx "
                       + std::to_string(orderArrayIdx) + ".");
        return nullptr;
      }
      case -3: {
        this->printErr("Retrieved enforced order array `"
                       + std::string(orderArray->GetName())
                       + "` has more than one component.");
        return nullptr;
      }
      case -2: {
        this->printErr("Enforced order array `"
                       + std::string(orderArray->GetName())
                       + "` is of incorrect type.");
        auto temp = vtkSmartPointer<ttkSimplexIdTypeArray>::New();
        this->printErr(" -> use `ttkArrayEditor` to convert data type to `"
                       + std::string(temp->GetDataTypeAsString()) + "`.");
        return nullptr;
      }
      default: {
        return checkForGlobalAndComputeOrderArray(
          inputData, scalarArray, scalarArrayIdx, getGlobalOrder, orderArray,
          triangulation, enforceOrderArrayIdx);
      }
    }
  }
 
  if(!scalarArray) {
    this->printErr("Unable to retrieve input scalar array for idx "
                   + std::to_string(scalarArrayIdx) + ".");
    return nullptr;
  } else if(isValidOrderArray(scalarArray) == 1) {
    this->printMsg("Retrieved scalar array `"
                     + std::string(scalarArray->GetName())
                     + "` is already an order array.",
                   ttk::debug::Priority::DETAIL);
    return scalarArray;
  }
 
  auto orderArray = inputData
                      ->GetAttributesAsFieldData(this->GetInputArrayAssociation(
                        scalarArrayIdx, inputData))
                      ->GetArray(this->GetOrderArrayName(scalarArray).data());
 
  switch(isValidOrderArray(orderArray)) {
    case -4: {
      ttk::Timer timer;
      printMsg(ttk::debug::Separator::L2);
      this->printWrn("No pre-existing order for array:");
      this->printWrn("  `" + std::string(scalarArray->GetName()) + "`.");
 
      this->printMsg("Initializing order array.", 0, 0, this->threadNumber_,
                     ttk::debug::LineMode::REPLACE);
      printMsg(ttk::debug::Separator::L2);
 
      orderArray
        = this->ComputeOrderArray(inputData, scalarArray, scalarArrayIdx,
                                  getGlobalOrder, orderArray, triangulation);
 
      std::string optionOn = "";
#ifdef TTK_ENABLE_MPI
      if(getGlobalOrder) {
        optionOn = "with GlobalOrder enabled ";
      }
      bool isGlobalOrder = getGlobalOrder || (!ttk::isRunningWithMPI());
      triangulation->setIsOrderArrayGlobal(
        ttkUtils::GetVoidPointer(scalarArray), isGlobalOrder);
#endif // TTK_ENABLE_MPI
      this->printMsg("Initializing order array.", 1, timer.getElapsedTime(),
                     this->threadNumber_);
 
      printMsg(ttk::debug::Separator::L2);
      this->printWrn("TIP: run `ttkArrayPreconditioning` first");
      this->printWrn(optionOn + "for improved performances :)");
      printMsg(ttk::debug::Separator::L2);
 
      return orderArray;
    }
 
    case -3: {
      this->printErr(
        "Retrieved order array `" + std::string(orderArray->GetName())
        + "` for scalar array `" + std::string(scalarArray->GetName())
        + "` has more than one component.");
      return nullptr;
    }
 
    case -2: {
      this->printErr(
        "Retrieved order array `" + std::string(orderArray->GetName())
        + "` for scalar array `" + std::string(scalarArray->GetName())
        + "` is of incorrect type.");
      auto temp = vtkSmartPointer<ttkSimplexIdTypeArray>::New();
      this->printErr(" -> use `ttkArrayEditor` to convert data type to `"
                     + std::string(temp->GetDataTypeAsString()) + "`.");
      return nullptr;
    }
 
    default: {
#ifdef TTK_ENABLE_MPI
      if(!ttk::isRunningWithMPI()) {
#endif
        this->printMsg("Retrieved order array `"
                         + std::string(orderArray->GetName())
                         + "` for scalar array `"
                         + std::string(scalarArray->GetName()) + "`.",
                       ttk::debug::Priority::DETAIL);
        return orderArray;
#ifdef TTK_ENABLE_MPI
      }
      return checkForGlobalAndComputeOrderArray(
        inputData, scalarArray, scalarArrayIdx, getGlobalOrder, orderArray,
        triangulation, enforceOrderArrayIdx);
#endif
    }
  }
}
 
ttk::SimplexId *
  ttkAlgorithm::GetIdentifierArrayPtr(const bool &enforceArrayIndex,
                                      const int &arrayIndex,
                                      const std::string &arrayName,
                                      vtkDataSet *const inputData,
                                      std::vector<ttk::SimplexId> &spareStorage,
                                      const int inputPort,
                                      const bool printErr) {
 
  // fetch data array
  const auto array = this->GetOptionalArray(
    enforceArrayIndex, arrayIndex, arrayName, inputData, inputPort);
  if(array == nullptr) {
    if(printErr) {
      this->printErr("Could not find the requested identifiers array");
    }
    return {};
  }
  if(array->GetNumberOfComponents() != 1) {
    if(printErr) {
      this->printErr("Identifiers field must have only one component!");
    }
    return {};
  }
 
#ifndef TTK_ENABLE_64BIT_IDS
  if(array->GetDataType() == VTK_ID_TYPE
     || array->GetDataType() == VTK_LONG_LONG) {
    this->printMsg(
      "Converting identifiers field from vtkIdType to SimplexId...");
    const auto nItems = array->GetNumberOfTuples();
 
    // fills the vector with the content of the data array converted to
    // ttk::SimplexId
    spareStorage.resize(nItems);
    for(vtkIdType i = 0; i < nItems; ++i) {
      spareStorage[i] = static_cast<ttk::SimplexId>(array->GetTuple1(i));
    }
 
    // return a pointer to the vector internal buffer
    return spareStorage.data();
  }
#else
  TTK_FORCE_USE(spareStorage);
#endif // TTK_ENABLE_64BIT_IDS
 
  // return a pointer to the data array internal buffer
  return static_cast<ttk::SimplexId *>(ttkUtils::GetVoidPointer(array));
}
 
template <class vtkDataType>
int prepOutput(vtkInformation *info, const std::string &className) {
  auto output = vtkDataObject::GetData(info);
  if(!output || !output->IsA(className.data())) {
    auto newOutput = vtkSmartPointer<vtkDataType>::New();
    info->Set(vtkDataObject::DATA_OBJECT(), newOutput);
  }
  return 1;
}
 
vtkDataSet *ttkAlgorithm::GetOutput() {
  return this->GetOutput(0);
}
 
vtkDataSet *ttkAlgorithm::GetOutput(int port) {
  return vtkDataSet::SafeDownCast(this->GetOutputDataObject(port));
}
 
void ttkAlgorithm::SetInputData(vtkDataSet *input) {
  this->SetInputData(0, input);
}
 
void ttkAlgorithm::SetInputData(int index, vtkDataSet *input) {
  this->SetInputDataInternal(index, input);
}
 
void ttkAlgorithm::AddInputData(vtkDataSet *input) {
  this->AddInputData(0, input);
}
 
void ttkAlgorithm::AddInputData(int index, vtkDataSet *input) {
  this->AddInputDataInternal(index, input);
}
 
int ttkAlgorithm::RequestDataObject(vtkInformation *ttkNotUsed(request),
                                    vtkInformationVector **inputVector,
                                    vtkInformationVector *outputVector) {
  // for each output
  for(int i = 0; i < this->GetNumberOfOutputPorts(); ++i) {
    auto outInfo = outputVector->GetInformationObject(i);
    if(!outInfo) {
      this->printErr("Unable to retrieve output vtkDataObject at port "
                     + std::to_string(i));
      return 0;
    }
 
    auto outputPortInfo = this->GetOutputPortInformation(i);
 
    // always request output type again for dynamic filter outputs
    if(!this->FillOutputPortInformation(i, outputPortInfo)) {
      this->printErr("Unable to fill output port information at port "
                     + std::to_string(i));
      return 0;
    }
 
    if(outputPortInfo->Has(ttkAlgorithm::SAME_DATA_TYPE_AS_INPUT_PORT())) {
      // Set output data type to input data type at specified port
      auto inPortIndex
        = outputPortInfo->Get(ttkAlgorithm::SAME_DATA_TYPE_AS_INPUT_PORT());
      if(inPortIndex < 0 || inPortIndex >= this->GetNumberOfInputPorts()) {
        this->printErr("Input port index " + std::to_string(inPortIndex)
                       + " specified by 'SAME_DATA_TYPE_AS_INPUT_PORT' key of "
                         "output port is out of range ("
                       + std::to_string(this->GetNumberOfInputPorts())
                       + " input ports).");
        return 0;
      }
      auto inInfo = inputVector[inPortIndex]->GetInformationObject(0);
      if(!inInfo) {
        this->printErr(
          "No information object at port " + std::to_string(inPortIndex)
          + " specified by 'SAME_DATA_TYPE_AS_INPUT_PORT' key of output port.");
        return 0;
      }
 
      auto input = vtkDataObject::GetData(inInfo);
      auto output = vtkDataObject::GetData(outInfo);
 
      if(!output || !output->IsA(input->GetClassName())) {
        auto newOutput
          = vtkSmartPointer<vtkDataObject>::Take(input->NewInstance());
        outputPortInfo->Set(
          vtkDataObject::DATA_TYPE_NAME(), input->GetClassName());
        outInfo->Set(vtkDataObject::DATA_OBJECT(), newOutput);
      }
    } else {
      // Explicitly create output by data type name
      if(!outputPortInfo->Has(vtkDataObject::DATA_TYPE_NAME())) {
        this->printErr("DATA_TYPE_NAME of output port " + std::to_string(i)
                       + " not specified");
        return 0;
      }
      std::string const outputType
        = outputPortInfo->Get(vtkDataObject::DATA_TYPE_NAME());
 
      if(outputType == "vtkUnstructuredGrid") {
        prepOutput<vtkUnstructuredGrid>(outInfo, outputType);
      } else if(outputType == "vtkPolyData") {
        prepOutput<vtkPolyData>(outInfo, outputType);
      } else if(outputType == "vtkMultiBlockDataSet") {
        prepOutput<vtkMultiBlockDataSet>(outInfo, outputType);
      } else if(outputType == "vtkTable") {
        prepOutput<vtkTable>(outInfo, outputType);
      } else if(outputType == "vtkImageData") {
        prepOutput<vtkImageData>(outInfo, outputType);
      } else {
        this->printErr("Unsupported data type for output[" + std::to_string(i)
                       + "]: " + outputType);
        return 0;
      }
    }
 
    this->printMsg(
      "Created '"
        + std::string(outputPortInfo->Get(vtkDataObject::DATA_TYPE_NAME()))
        + "' at output port " + std::to_string(i),
      ttk::debug::Priority::VERBOSE);
  }
 
  return 1;
}
 
#ifdef TTK_ENABLE_MPI
 
int ttkAlgorithm::updateMPICommunicator(vtkDataSet *input) {
  if(input == nullptr) {
    return 0;
  }
  int isEmpty
    = input->GetNumberOfCells() == 0 || input->GetNumberOfPoints() == 0;
  int oldSize = ttk::MPIsize_;
  int oldRank = ttk::MPIrank_;
  MPI_Comm_split(MPI_COMM_WORLD, isEmpty, 0, &ttk::MPIcomm_);
  MPI_Comm_rank(ttk::MPIcomm_, &ttk::MPIrank_);
  MPI_Comm_size(ttk::MPIcomm_, &ttk::MPIsize_);
  if(oldSize != ttk::MPIsize_) {
    std::vector<int> newToOldRanks(ttk::MPIsize_);
    MPI_Allgather(&oldRank, 1, MPI_INTEGER, newToOldRanks.data(), 1,
                  MPI_INTEGER, ttk::MPIcomm_);
    std::map<int, int> oldToNewRanks;
    for(int i = 0; i < ttk::MPIsize_; i++) {
      oldToNewRanks[newToOldRanks[i]] = i;
    }
    int *vertexRankArray
      = ttkUtils::GetPointer<int>(input->GetPointData()->GetArray("RankArray"));
    if(vertexRankArray != nullptr) {
      for(int i = 0; i < input->GetNumberOfPoints(); i++) {
        vertexRankArray[i] = oldToNewRanks[vertexRankArray[i]];
      }
    }
    int *cellRankArray
      = ttkUtils::GetPointer<int>(input->GetCellData()->GetArray("RankArray"));
    if(cellRankArray != nullptr) {
      for(int i = 0; i < input->GetNumberOfCells(); i++) {
        cellRankArray[i] = oldToNewRanks[cellRankArray[i]];
      }
    }
  }
  this->setDebugMsgPrefix(debugMsgNamePrefix_);
  return isEmpty;
}
 
bool ttkAlgorithm::checkGlobalIdValidity(ttk::LongSimplexId *globalIds,
                                         ttk::SimplexId simplexNumber,
                                         unsigned char *ghost,
                                         int *rankArray) {
  ttk::SimplexId ghostNumber = 0;
  if(rankArray != nullptr) {
#ifdef TTK_ENABLE_OPENMP
#pragma omp parallel for reduction(+ : ghostNumber)
#endif // TTK_ENABLE_OPENMP
    for(ttk::SimplexId i = 0; i < simplexNumber; i++) {
      if(rankArray[i] != ttk::MPIrank_) {
        ghostNumber++;
      }
    }
  } else {
    if(ghost != nullptr) {
#ifdef TTK_ENABLE_OPENMP
#pragma omp parallel for reduction(+ : ghostNumber)
#endif // TTK_ENABLE_OPENMP
      for(ttk::SimplexId i = 0; i < simplexNumber; i++) {
        if(ghost[i] == 1) {
          ghostNumber++;
        }
      }
    }
  }
 
  ttk::SimplexId realSimplexNumber = simplexNumber - ghostNumber;
  auto minmax = std::minmax_element(globalIds, globalIds + simplexNumber);
  ttk::LongSimplexId min = globalIds[minmax.first - globalIds];
  ttk::LongSimplexId max = globalIds[minmax.second - globalIds];
  ttk::SimplexId globalSimplexNumber;
  ttk::LongSimplexId globalMin;
  ttk::LongSimplexId globalMax;
  MPI_Allreduce(&realSimplexNumber, &globalSimplexNumber, 1,
                ttk::getMPIType(realSimplexNumber), MPI_SUM, ttk::MPIcomm_);
  MPI_Allreduce(
    &min, &globalMin, 1, ttk::getMPIType(min), MPI_MIN, ttk::MPIcomm_);
  MPI_Allreduce(
    &max, &globalMax, 1, ttk::getMPIType(max), MPI_MAX, ttk::MPIcomm_);
 
  return (globalSimplexNumber == globalMax + 1 && globalMin == 0);
};
 
int ttkAlgorithm::GenerateGlobalIds(
  vtkDataSet *input,
  std::unordered_map<ttk::SimplexId, ttk::SimplexId> &vertGtoL,
  std::vector<int> &neighborRanks,
  std::map<int, int> &neighborsToId) {
 
  ttk::Identifiers identifiers;
 
  vtkNew<vtkIdTypeArray> vtkVertexIdentifiers{};
 
  vtkNew<vtkIdTypeArray> vtkCellIdentifiers{};
  vtkVertexIdentifiers->SetName("GlobalPointIds");
  vtkVertexIdentifiers->SetNumberOfComponents(1);
  vtkVertexIdentifiers->SetNumberOfTuples(input->GetNumberOfPoints());
  vtkVertexIdentifiers->Fill(-1);
 
  identifiers.setVertexIdentifiers(
    ttkUtils::GetPointer<ttk::LongSimplexId>(vtkVertexIdentifiers));
 
  vtkCellIdentifiers->SetName("GlobalCellIds");
  vtkCellIdentifiers->SetNumberOfComponents(1);
  vtkCellIdentifiers->SetNumberOfTuples(input->GetNumberOfCells());
  vtkCellIdentifiers->Fill(-1);
 
  identifiers.setCellIdentifiers(
    ttkUtils::GetPointer<ttk::LongSimplexId>(vtkCellIdentifiers));
 
  int vertexNumber = input->GetNumberOfPoints();
  identifiers.setVertexNumber(vertexNumber);
  int cellNumber = input->GetNumberOfCells();
  identifiers.setCellNumber(cellNumber);
  int status = 0;
 
  double *boundingBox = input->GetBounds();
  identifiers.setBounds(boundingBox);
  identifiers.initializeNeighbors(boundingBox, neighborRanks, neighborsToId);
  if(ttk::isRunningWithMPI()) {
    switch(input->GetDataObjectType()) {
      case VTK_UNSTRUCTURED_GRID:
      case VTK_POLY_DATA: {
 
        identifiers.setOutdatedGlobalPointIds(
          ttkUtils::GetPointer<ttk::LongSimplexId>(
            input->GetPointData()->GetGlobalIds()));
        identifiers.setOutdatedGlobalCellIds(
          ttkUtils::GetPointer<ttk::LongSimplexId>(
            input->GetCellData()->GetGlobalIds()));
        identifiers.setVertexRankArray(ttkUtils::GetPointer<int>(
          input->GetPointData()->GetArray("RankArray")));
        int *cellRankArray = ttkUtils::GetPointer<int>(
          input->GetCellData()->GetArray("RankArray"));
        identifiers.setCellRankArray(cellRankArray);
        identifiers.setVertGhost(ttkUtils::GetPointer<unsigned char>(
          input->GetPointData()->GetArray("vtkGhostType")));
        unsigned char *cellGhost = ttkUtils::GetPointer<unsigned char>(
          input->GetCellData()->GetArray("vtkGhostType"));
        identifiers.setCellGhost(cellGhost);
        vtkPointSet *pointSet = vtkPointSet::SafeDownCast(input);
        identifiers.setPointSet(static_cast<float *>(
          ttkUtils::GetVoidPointer(pointSet->GetPoints())));
        vtkCellArray *cells = nullptr;
        switch(input->GetDataObjectType()) {
          case VTK_UNSTRUCTURED_GRID: {
            auto dataSetAsUG = vtkUnstructuredGrid::SafeDownCast(input);
            cells = dataSetAsUG->GetCells();
            break;
          }
          case VTK_POLY_DATA: {
            auto dataSetAsPD = vtkPolyData::SafeDownCast(input);
            cells
              = dataSetAsPD->GetNumberOfPolys() > 0   ? dataSetAsPD->GetPolys()
                : dataSetAsPD->GetNumberOfLines() > 0 ? dataSetAsPD->GetLines()
                                                      : dataSetAsPD->GetVerts();
            break;
          }
          default: {
            this->printErr("Unable to get cells for `"
                           + std::string(input->GetClassName()) + "`");
          }
        }
        if(cells == nullptr) {
          return 0;
        }
        if(!cells->IsStorage64Bit()) {
          if(cells->CanConvertTo64BitStorage()) {
            this->printWrn("Converting the cell array to 64-bit storage");
            bool success = cells->ConvertTo64BitStorage();
            if(!success) {
              this->printErr(
                "Error converting the provided cell array to 64-bit storage");
              return -1;
            }
          } else {
            this->printErr(
              "Cannot convert the provided cell array to 64-bit storage");
            return -1;
          }
        }
 
        identifiers.setConnectivity(ttkUtils::GetPointer<ttk::LongSimplexId>(
          cells->GetConnectivityArray()));
 
        std::vector<std::vector<ttk::SimplexId>> pointsToCells(vertexNumber);
        vtkIdList *cellList = vtkIdList::New();
        if(cellRankArray != nullptr) {
          for(ttk::SimplexId i = 0; i < vertexNumber; i++) {
            input->GetPointCells(i, cellList);
            for(int j = 0; j < cellList->GetNumberOfIds(); j++) {
              if(cellRankArray[cellList->GetId(j)] == ttk::MPIrank_) {
                pointsToCells[i].push_back(cellList->GetId(j));
              }
            }
          }
        } else {
          for(ttk::SimplexId i = 0; i < vertexNumber; i++) {
            input->GetPointCells(i, cellList);
            for(int j = 0; j < cellList->GetNumberOfIds(); j++) {
              if(cellGhost[cellList->GetId(j)] == 0) {
                pointsToCells[i].push_back(cellList->GetId(j));
              }
            }
          }
        }
        identifiers.setPointsToCells(pointsToCells);
 
        identifiers.initializeMPITypes();
        identifiers.setVertGtoL(&vertGtoL);
        vtkIdList *pointCell = vtkIdList::New();
        input->GetCellPoints(0, pointCell);
        int nbPoints = pointCell->GetNumberOfIds();
        identifiers.setDomainDimension(nbPoints - 1);
        identifiers.buildKDTree();
        status = identifiers.executePolyData();
        break;
      }
      case VTK_IMAGE_DATA: {
        vtkImageData *data = vtkImageData::SafeDownCast(input);
        identifiers.setDims(data->GetDimensions());
        identifiers.setSpacing(data->GetSpacing());
        status = identifiers.executeImageData();
        break;
      }
      default: {
        this->printErr("Unable to triangulate `"
                       + std::string(input->GetClassName()) + "`");
      }
    }
  } else {
    status = identifiers.executeSequential();
  }
 
  if(status < 1) {
    printErr("Global identifier generation failed");
    return -1;
  }
 
  // Add VTK objects to the data set
 
  input->GetPointData()->SetGlobalIds(vtkVertexIdentifiers);
  input->GetCellData()->SetGlobalIds(vtkCellIdentifiers);
  return 0;
}
 
void ttkAlgorithm::MPIGhostPipelinePreconditioning(vtkDataSet *input) {
 
  vtkNew<vtkGhostCellsGenerator> generator;
  if(ttk::isRunningWithMPI()
     && (!input->HasAnyGhostCells()
         && ((input->GetPointData()->GetArray("RankArray") == nullptr)
             || (input->GetCellData()->GetArray("RankArray") == nullptr)))) {
    generator->SetInputData(input);
    generator->BuildIfRequiredOff();
    generator->SetNumberOfGhostLayers(1);
    generator->Update();
    input->ShallowCopy(generator->GetOutputDataObject(0));
    input->GetPointData()->AddArray(
      generator->GetOutputDataObject(0)->GetGhostArray(0));
    input->GetCellData()->AddArray(
      generator->GetOutputDataObject(0)->GetGhostArray(1));
  }
}
 
void ttkAlgorithm::MPIPipelinePreconditioning(
  vtkDataSet *input,
  std::vector<int> &neighbors,
  std::map<int, int> &neighToId,
  ttk::Triangulation *triangulation) {
 
  ttk::SimplexId vertexNumber = input->GetNumberOfPoints();
  ttk::SimplexId cellNumber = input->GetNumberOfCells();
 
  if((input->GetDataObjectType() == VTK_POLY_DATA
      || input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)) {
 
    if((ttk::hasInitializedMPI()) && (ttk::isRunningWithMPI())) {
      printMsg(ttk::debug::Separator::L2);
      printWrn("The distribution by VTK of Unstructured");
      printWrn("Grids and Poly Data has been reported");
      printWrn("to be affected by bugs (at least up");
      printWrn("to ParaView 5.10.1).");
 
      if((input->GetCellData()->GetGlobalIds() == nullptr)
         || (input->GetPointData()->GetGlobalIds() == nullptr)) {
 
        printWrn("=> Global identifiers may be incorrect.");
      }
      if((input->GetPointData()->GetArray("RankArray") == nullptr)
         || (input->GetCellData()->GetArray("RankArray") == nullptr)) {
        printWrn("=> Rank arrays may be incorrect.");
      }
      printMsg(ttk::debug::Separator::L2);
    }
  }
 
  // Get the neighbor ranks
  std::vector<int> &neighborRanks{
    triangulation != nullptr ? triangulation->getNeighborRanks() : neighbors};
  std::map<int, int> &neighborsToId{
    triangulation != nullptr ? triangulation->getNeighborsToId() : neighToId};
 
  double *boundingBox = input->GetBounds();
  if(triangulation != nullptr) {
    triangulation->createMetaGrid(boundingBox);
  }
 
  if(neighborRanks.empty()) {
    ttk::preconditionNeighborsUsingBoundingBox(
      boundingBox, neighborRanks, neighborsToId);
  }
 
  // Checks if global ids are valid
  ttk::LongSimplexId *globalPointIds = ttkUtils::GetPointer<ttk::LongSimplexId>(
    input->GetPointData()->GetGlobalIds());
  ttk::LongSimplexId *globalCellIds = ttkUtils::GetPointer<ttk::LongSimplexId>(
    input->GetCellData()->GetGlobalIds());
 
  bool pointValidity{false};
  bool cellValidity{false};
  if((triangulation != nullptr
      && (triangulation->getType() == ttk::Triangulation::Type::EXPLICIT
          || triangulation->getType() == ttk::Triangulation::Type::COMPACT))
     || triangulation == nullptr) {
    if(globalPointIds != nullptr) {
      unsigned char *ghostPoints = ttkUtils::GetPointer<unsigned char>(
        input->GetPointData()->GetArray("vtkGhostType"));
      int *vertexRankArray = ttkUtils::GetPointer<int>(
        input->GetPointData()->GetArray("RankArray"));
      pointValidity = checkGlobalIdValidity(
        globalPointIds, vertexNumber, ghostPoints, vertexRankArray);
    }
    if(pointValidity && globalCellIds != nullptr) {
 
      unsigned char *ghostCells = ttkUtils::GetPointer<unsigned char>(
        input->GetCellData()->GetArray("vtkGhostType"));
      int *cellRankArray = ttkUtils::GetPointer<int>(
        input->GetCellData()->GetArray("RankArray"));
      cellValidity = checkGlobalIdValidity(
        globalCellIds, cellNumber, ghostCells, cellRankArray);
    }
  } else {
    pointValidity = true;
    cellValidity = true;
  }
 
  // If the global ids are not valid, they are computed again
  if(!pointValidity || !cellValidity) {
    if(triangulation != nullptr) {
      if(triangulation->getType() == ttk::Triangulation::Type::EXPLICIT
         || triangulation->getType() == ttk::Triangulation::Type::COMPACT) {
        this->GenerateGlobalIds(input, triangulation->getVertexGlobalIdMap(),
                                neighborRanks, neighborsToId);
      }
    } else {
      std::unordered_map<ttk::SimplexId, ttk::SimplexId> vertGtoL{};
      this->GenerateGlobalIds(input, vertGtoL, neighborRanks, neighborsToId);
    }
  }
}
 
void ttkAlgorithm::MPITriangulationPreconditioning(
  ttk::Triangulation *triangulation, vtkDataSet *input) {
 
  const auto pd{input->GetPointData()};
  if(pd == nullptr) {
    triangulation->printWrn("No point data on input object");
  } else {
    // provide "GlobalPointIds" & "vtkGhostType" point data arrays
    // to the triangulation
    triangulation->setVertsGlobalIds(
      ttkUtils::GetPointer<ttk::LongSimplexId>(pd->GetGlobalIds()));
    triangulation->setVertexGhostArray(
      ttkUtils::GetPointer<unsigned char>(pd->GetArray("vtkGhostType")));
    int *vertexRankArray = ttkUtils::GetPointer<int>(pd->GetArray("RankArray"));
    if(vertexRankArray != nullptr) {
      triangulation->setVertexRankArray(vertexRankArray);
    }
    triangulation->preconditionDistributedVertices();
  }
 
  const auto cd{input->GetCellData()};
  if(cd == nullptr) {
    triangulation->printWrn("No cell data on input object");
  } else {
    // provide "GlobalCellIds" & "vtkGhostType" cell data arrays to
    // the triangulation
    triangulation->setCellsGlobalIds(
      ttkUtils::GetPointer<ttk::LongSimplexId>(cd->GetGlobalIds()));
    triangulation->setCellGhostArray(
      ttkUtils::GetPointer<unsigned char>(cd->GetArray("vtkGhostType")));
    int *cellRankArray = ttkUtils::GetPointer<int>(cd->GetArray("RankArray"));
    if(cellRankArray != nullptr) {
      triangulation->setCellRankArray(cellRankArray);
    }
    triangulation->preconditionDistributedCells();
  }
}
 
#endif // TTK_ENABLE_MPI
 
//==============================================================================
int ttkAlgorithm::ProcessRequest(vtkInformation *request,
                                 vtkInformationVector **inputVector,
                                 vtkInformationVector *outputVector) {
  // 1. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA_OBJECT())) {
    this->printMsg(
      "Processing REQUEST_DATA_OBJECT", ttk::debug::Priority::VERBOSE);
    return this->RequestDataObject(request, inputVector, outputVector);
  }
 
  // 2. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_INFORMATION())) {
    this->printMsg(
      "Processing REQUEST_INFORMATION", ttk::debug::Priority::VERBOSE);
    return this->RequestInformation(request, inputVector, outputVector);
  }
 
  // 3. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_UPDATE_TIME())) {
    this->printMsg(
      "Processing REQUEST_UPDATE_TIME", ttk::debug::Priority::VERBOSE);
    return this->RequestUpdateTime(request, inputVector, outputVector);
  }
 
  // 4. Pass
  if(request->Has(
       vtkCompositeDataPipeline::REQUEST_TIME_DEPENDENT_INFORMATION())) {
    this->printMsg("Processing REQUEST_TIME_DEPENDENT_INFORMATION",
                   ttk::debug::Priority::VERBOSE);
    return this->RequestUpdateTimeDependentInformation(
      request, inputVector, outputVector);
  }
 
  // 5. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_UPDATE_EXTENT())) {
    this->printMsg(
      "Processing REQUEST_UPDATE_EXTENT", ttk::debug::Priority::VERBOSE);
    return this->RequestUpdateExtent(request, inputVector, outputVector);
  }
 
  // 6. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA_NOT_GENERATED())) {
    this->printMsg(
      "Processing REQUEST_DATA_NOT_GENERATED", ttk::debug::Priority::VERBOSE);
    return this->RequestDataNotGenerated(request, inputVector, outputVector);
  }
 
  // 7. Pass
  if(request->Has(vtkCompositeDataPipeline::REQUEST_DATA())) {
    this->printMsg("Processing REQUEST_DATA", ttk::debug::Priority::VERBOSE);
    this->printMsg(ttk::debug::Separator::L0);
#ifdef TTK_ENABLE_MPI
    if(ttk::hasInitializedMPI() && inputVector != nullptr) {
      if(this->updateMPICommunicator(vtkDataSet::GetData(inputVector[0], 0))) {
        return 1;
      };
    }
#endif // TTK_ENABLE_MPI
    return this->RequestData(request, inputVector, outputVector);
  }
 
  this->printErr("Unsupported pipeline pass:");
  request->Print(cout);
 
  return 0;
}