text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    vtkDistributedDataFilter.cxx
  
    Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
    All rights reserved.
    See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
  
       This software is distributed WITHOUT ANY WARRANTY; without even
       the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
       PURPOSE.  See the above copyright notice for more information.
  
  =========================================================================*/
  /*----------------------------------------------------------------------------
   Copyright (c) Sandia Corporation
   See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
  ----------------------------------------------------------------------------*/
  
  #include "vtkDistributedDataFilter.h"
  
  #include "vtkBox.h"
  #include "vtkBoxClipDataSet.h"
  #include "vtkBSPCuts.h"
  #include "vtkCellArray.h"
  #include "vtkCellData.h"
  #include "vtkCharArray.h"
  #include "vtkClipDataSet.h"
  #include "vtkCompositeDataIterator.h"
  #include "vtkDataObjectTypes.h"
  #include "vtkDataSetReader.h"
  #include "vtkDataSetWriter.h"
  #include "vtkExtractCells.h"
  #include "vtkExtractUserDefinedPiece.h"
  #include "vtkFloatArray.h"
  #include "vtkIdList.h"
  #include "vtkIdTypeArray.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkIntArray.h"
  #include "vtkMergeCells.h"
  #include "vtkMultiBlockDataSet.h"
  #include "vtkMultiProcessController.h"
  #include "vtkObjectFactory.h"
  #include "vtkPKdTree.h"
  #include "vtkPlane.h"
  #include "vtkPointData.h"
  #include "vtkPointLocator.h"
  #include "vtkSmartPointer.h"
  #include "vtkSocketController.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  #include "vtkToolkits.h"
  #include "vtkUnsignedCharArray.h"
  #include "vtkUnstructuredGrid.h"
  
  #include "vtkMPIController.h"
  
  #include <vector>
  
  
  vtkStandardNewMacro(vtkDistributedDataFilter)
  
  #define TEMP_ELEMENT_ID_NAME      "___D3___GlobalCellIds"
  #define TEMP_INSIDE_BOX_FLAG      "___D3___WHERE"
  #define TEMP_NODE_ID_NAME         "___D3___GlobalNodeIds"
  
  #include <set>
  #include <map>
  #include <algorithm>
  
  class vtkDistributedDataFilterSTLCloak
  {
  public:
    std::map<int, int> IntMap;
    std::multimap<int, int> IntMultiMap;
  };
  
  namespace
  {
  void convertGhostLevelsToBitFields(vtkDataSetAttributes* dsa, unsigned int bit)
  {
    vtkDataArray *da = dsa->GetArray(vtkDataSetAttributes::GhostArrayName());
    vtkUnsignedCharArray *uca = vtkArrayDownCast<vtkUnsignedCharArray>(da);
    unsigned char *ghosts = uca->GetPointer(0);
    for (vtkIdType i=0; i < da->GetNumberOfTuples(); ++i)
    {
      if (ghosts[i] > 0)
      {
        ghosts[i] = bit;
      }
    }
  }
  };
  
  class vtkDistributedDataFilter::vtkInternals
  {
  public:
    std::vector<int> UserRegionAssignments;
  };
  
  //----------------------------------------------------------------------------
  vtkDistributedDataFilter::vtkDistributedDataFilter()
  {
    this->Kdtree = nullptr;
  
    this->Controller = nullptr;
    this->SetController(vtkMultiProcessController::GetGlobalController());
  
    this->Target = nullptr;
    this->Source = nullptr;
  
    this->NumConvexSubRegions = 0;
    this->ConvexSubRegionBounds = nullptr;
  
    this->GhostLevel = 0;
  
    this->RetainKdtree = 1;
    this->IncludeAllIntersectingCells = 0;
    this->ClipCells = 0;
  
    this->Timing = 0;
  
    this->UseMinimalMemory = 0;
  
    this->UserCuts = nullptr;
    this->Internals = new vtkDistributedDataFilter::vtkInternals();
  }
  
  //----------------------------------------------------------------------------
  vtkDistributedDataFilter::~vtkDistributedDataFilter()
  {
    if (this->Kdtree)
    {
      this->Kdtree->Delete();
      this->Kdtree = nullptr;
    }
  
    this->SetController(nullptr);
  
    delete [] this->Target;
    this->Target= nullptr;
  
    delete [] this->Source;
    this->Source= nullptr;
  
    delete [] this->ConvexSubRegionBounds;
    this->ConvexSubRegionBounds = nullptr;
  
    if (this->UserCuts)
    {
      this->UserCuts->Delete();
      this->UserCuts = nullptr;
    }
    delete this->Internals;
    this->Internals = nullptr;
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetCuts(vtkBSPCuts* cuts)
  {
    if (cuts == this->UserCuts)
    {
      return;
    }
    if (this->UserCuts)
    {
      this->UserCuts->Delete();
      this->UserCuts = nullptr;
    }
    if (cuts)
    {
      cuts->Register(this);
      this->UserCuts = cuts;
    }
    // Delete the Kdtree so that it is regenerated next time.
    if (this->Kdtree)
    {
      this->Kdtree->SetCuts(cuts);
    }
    this->Modified();
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetUserRegionAssignments(
    const int *map, int numRegions)
  {
    std::vector<int> copy(this->Internals->UserRegionAssignments);
    this->Internals->UserRegionAssignments.resize(numRegions);
    for (int cc=0; cc < numRegions; cc++)
    {
      this->Internals->UserRegionAssignments[cc] = map[cc];
    }
    if (copy != this->Internals->UserRegionAssignments)
    {
      this->Modified();
    }
  }
  
  //----------------------------------------------------------------------------
  vtkIdTypeArray *vtkDistributedDataFilter::GetGlobalElementIdArray(vtkDataSet *set)
  {
    vtkDataArray *da = set->GetCellData()->GetGlobalIds();
    return vtkArrayDownCast<vtkIdTypeArray>(da);
  }
  
  //----------------------------------------------------------------------------
  vtkIdType *vtkDistributedDataFilter::GetGlobalElementIds(vtkDataSet *set)
  {
    vtkIdTypeArray *ia = GetGlobalElementIdArray(set);
    if (!ia)
    {
      return nullptr;
    }
  
    return ia->GetPointer(0);
  }
  
  //----------------------------------------------------------------------------
  vtkIdTypeArray *vtkDistributedDataFilter::GetGlobalNodeIdArray(vtkDataSet *set)
  {
    vtkDataArray *da = set->GetPointData()->GetGlobalIds();
    return vtkArrayDownCast<vtkIdTypeArray>(da);
  }
  
  //----------------------------------------------------------------------------
  vtkIdType *vtkDistributedDataFilter::GetGlobalNodeIds(vtkDataSet *set)
  {
    vtkIdTypeArray *ia = this->GetGlobalNodeIdArray(set);
  
    if (!ia)
    {
      return nullptr;
    }
  
    return ia->GetPointer(0);
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetController(vtkMultiProcessController *c)
  {
    if (this->Kdtree)
    {
      this->Kdtree->SetController(c);
    }
  
    if ((c == nullptr) || (c->GetNumberOfProcesses() == 0))
    {
      this->NumProcesses = 1;
      this->MyId = 0;
    }
  
    if (this->Controller == c)
    {
      return;
    }
  
    this->Modified();
  
    if (this->Controller != nullptr)
    {
      this->Controller->UnRegister(this);
      this->Controller = nullptr;
    }
  
    if (c == nullptr)
    {
      return;
    }
  
    this->Controller = c;
  
    c->Register(this);
    this->NumProcesses = c->GetNumberOfProcesses();
    this->MyId    = c->GetLocalProcessId();
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetBoundaryMode(int mode)
  {
    switch (mode)
    {
      case vtkDistributedDataFilter::ASSIGN_TO_ONE_REGION:
        this->AssignBoundaryCellsToOneRegionOn();
        break;
      case vtkDistributedDataFilter::ASSIGN_TO_ALL_INTERSECTING_REGIONS:
        this->AssignBoundaryCellsToAllIntersectingRegionsOn();
        break;
      case vtkDistributedDataFilter::SPLIT_BOUNDARY_CELLS:
        this->DivideBoundaryCellsOn();
        break;
    }
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::GetBoundaryMode()
  {
    if (!this->IncludeAllIntersectingCells && !this->ClipCells)
    {
      return vtkDistributedDataFilter::ASSIGN_TO_ONE_REGION;
    }
    if (this->IncludeAllIntersectingCells && !this->ClipCells)
    {
      return vtkDistributedDataFilter::ASSIGN_TO_ALL_INTERSECTING_REGIONS;
    }
    if (this->IncludeAllIntersectingCells && this->ClipCells)
    {
      return vtkDistributedDataFilter::SPLIT_BOUNDARY_CELLS;
    }
  
    return -1;
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::AssignBoundaryCellsToOneRegionOn()
  {
    this->SetAssignBoundaryCellsToOneRegion(1);
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::AssignBoundaryCellsToOneRegionOff()
  {
    this->SetAssignBoundaryCellsToOneRegion(0);
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetAssignBoundaryCellsToOneRegion(int val)
  {
    if (val)
    {
      this->IncludeAllIntersectingCells = 0;
      this->ClipCells = 0;
    }
  }
  
  //----------------------------------------------------------------------------
  void
  vtkDistributedDataFilter::AssignBoundaryCellsToAllIntersectingRegionsOn()
  {
    this->SetAssignBoundaryCellsToAllIntersectingRegions(1);
  }
  
  //----------------------------------------------------------------------------
  void
  vtkDistributedDataFilter::AssignBoundaryCellsToAllIntersectingRegionsOff()
  {
    this->SetAssignBoundaryCellsToAllIntersectingRegions(0);
  }
  
  //----------------------------------------------------------------------------
  void
  vtkDistributedDataFilter::SetAssignBoundaryCellsToAllIntersectingRegions(int val)
  {
    if (val)
    {
      this->IncludeAllIntersectingCells = 1;
      this->ClipCells = 0;
    }
  }
  void vtkDistributedDataFilter::DivideBoundaryCellsOn()
  {
    this->SetDivideBoundaryCells(1);
  }
  void vtkDistributedDataFilter::DivideBoundaryCellsOff()
  {
    this->SetDivideBoundaryCells(0);
  }
  void vtkDistributedDataFilter::SetDivideBoundaryCells(int val)
  {
    if (val)
    {
      this->IncludeAllIntersectingCells = 1;
      this->ClipCells = 1;
    }
  }
  
  
  //============================================================================
  // Execute
  
  //-------------------------------------------------------------------------
  
  int vtkDistributedDataFilter::RequestUpdateExtent(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    // get the info objects
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    int piece, numPieces, ghostLevels;
  
    // We require preceding filters to refrain from creating ghost cells.
  
    piece = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
    numPieces = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
    ghostLevels = 0;
  
    inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), piece);
    inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
                numPieces);
    inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
                ghostLevels);
    inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);
  
    return 1;
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::RequestInformation(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    // get the info objects
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
                 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
                 6);
  
    return 1;
  }
  
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    // get the info objects
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    this->GhostLevel = outInfo->Get(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
  
    // get the input and output
    vtkDataSet *inputDS = vtkDataSet::GetData(inputVector[0], 0);
    vtkUnstructuredGrid *outputUG = vtkUnstructuredGrid::GetData(outInfo);
    if (inputDS && outputUG)
    {
      return this->RequestDataInternal(inputDS, outputUG);
    }
  
    vtkCompositeDataSet* inputCD =
      vtkCompositeDataSet::GetData(inputVector[0], 0);
    vtkMultiBlockDataSet* outputMB =
      vtkMultiBlockDataSet::GetData(outputVector, 0);
    if (!inputCD || !outputMB)
    {
      vtkErrorMacro("Input must either be a composite dataset or a vtkDataSet.");
      return 0;
    }
  
    outputMB->CopyStructure(inputCD);
  
    vtkSmartPointer<vtkCompositeDataIterator> iter;
    iter.TakeReference(inputCD->NewIterator());
    // We want to traverse over empty nodes as well. This ensures that this
    // algorithm will work correctly in parallel.
    iter->SkipEmptyNodesOff();
  
    // Collect information about datatypes all the processes have at all the leaf
    // nodes. Ideally all processes will either have the same type or an empty
    // dataset. This assumes that all processes have the same composite structure.
    std::vector<int> leafTypes;
    for (iter->InitTraversal(); !iter->IsDoneWithTraversal();
      iter->GoToNextItem())
    {
      vtkDataObject* dObj = iter->GetCurrentDataObject();
      if (dObj)
      {
        leafTypes.push_back(dObj->GetDataObjectType());
      }
      else
      {
        leafTypes.push_back(-1);
      }
    }
    unsigned int numLeaves = static_cast<unsigned int>(leafTypes.size());
  
    int myId = this->Controller->GetLocalProcessId();
    int numProcs = this->Controller->GetNumberOfProcesses();
    if (numProcs > 1 && numLeaves > 0)
    {
      if (myId == 0)
      {
        for (int cc=1; cc < numProcs; cc++)
        {
          std::vector<int> receivedTypes;
          receivedTypes.resize(numLeaves, -1);
          if (!this->Controller->Receive(&receivedTypes[0],
              numLeaves, cc, 1020202))
          {
            vtkErrorMacro("Communication error.");
            return 0;
          }
          for (unsigned int kk=0; kk < numLeaves; kk++)
          {
            if (leafTypes[kk] == -1)
            {
              leafTypes[kk] = receivedTypes[kk];
            }
            if (receivedTypes[kk] != -1 && leafTypes[kk] != -1 &&
              receivedTypes[kk] != leafTypes[kk])
            {
              vtkWarningMacro("Data type mismatch on processes.");
            }
          }
        }
        for (int kk=1; kk < numProcs; kk++)
        {
          this->Controller->Send(&leafTypes[0], numLeaves, kk, 1020203);
        }
      }
      else
      {
        this->Controller->Send(&leafTypes[0], numLeaves, 0, 1020202);
        this->Controller->Receive(&leafTypes[0], numLeaves, 0, 1020203);
      }
    }
  
    unsigned int cc=0;
    for (iter->InitTraversal(); !iter->IsDoneWithTraversal();
      iter->GoToNextItem(), cc++)
    {
      vtkSmartPointer<vtkDataSet> ds;
      ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
      if (ds.GetPointer() == nullptr)
      {
        if (leafTypes[cc] == -1)
        {
          // This is an empty block on all processes, just skip it.
          continue;
        }
  
        ds.TakeReference(vtkDataSet::SafeDownCast(
          vtkDataObjectTypes::NewDataObject(leafTypes[cc])));
      }
      vtkSmartPointer<vtkUnstructuredGrid> ug =
        vtkSmartPointer<vtkUnstructuredGrid>::New();
      if (!this->RequestDataInternal(ds, ug))
      {
        return 0;
      }
      if (ug->GetNumberOfPoints() >0)
      {
        outputMB->SetDataSet(iter, ug);
      }
    }
    return 1;
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::RequestDataInternal(vtkDataSet* input,
    vtkUnstructuredGrid* output)
  {
    this->NextProgressStep = 0;
    int progressSteps = 5 + this->GhostLevel;
    if (this->ClipCells)
    {
      progressSteps++;
    }
  
    this->ProgressIncrement = 1.0 / (double)progressSteps;
  
    this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
    this->SetProgressText("Begin data redistribution");
  
    if (this->NumProcesses == 1)
    {
      this->SingleProcessExecute(input, output);
      this->UpdateProgress(1.0);
      return 1;
    }
  
    // This method requires an MPI controller.
  
    int aok = 0;
  
    if (vtkMPIController::SafeDownCast(this->Controller))
    {
      aok = 1;
    }
  
    if (!aok)
    {
      vtkErrorMacro(<< "vtkDistributedDataFilter multiprocess requires MPI");
      return 1;
    }
  
    // Stage (0) - If any processes have 0 cell input data sets, then
    //   spread the input data sets around (quickly) before formal
    //   redistribution.
  
    vtkDataSet *splitInput = this->TestFixTooFewInputFiles(input);
  
    if (splitInput == nullptr)
    {
      return 1;    // Fewer cells than processes - can't divide input
    }
  
    this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
    this->SetProgressText("Compute spatial partitioning");
  
    // Stage (1) - use vtkPKdTree to...
    //   Create a load balanced spatial decomposition in parallel.
    //   Create a table assigning regions to processes.
    //
    // Note k-d tree will only be re-built if input or parameters
    // have changed on any of the processing nodes.
  
    int fail = this->PartitionDataAndAssignToProcesses(splitInput);
  
    if (fail)
    {
      if (splitInput != input)
      {
        splitInput->Delete();
      }
      vtkErrorMacro(<< "vtkDistributedDataFilter::Execute k-d tree failure");
      return 1;
    }
  
    this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
    this->SetProgressText("Compute global data array bounds");
  
    // Let the vtkPKdTree class compile global bounds for all
    // data arrays.  These can be accessed by D3 user by getting
    // a handle to the vtkPKdTree object and querying it.
  
    this->Kdtree->CreateGlobalDataArrayBounds();
  
    this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
    this->SetProgressText("Redistribute data");
  
    // Stage (2) - Redistribute data, so that each process gets a ugrid
    //   containing the cells in it's assigned spatial regions.  (Note
    //   that a side effect of merging the grids received from different
    //   processes is that the final grid has no duplicate points.)
    //
    // This call will delete splitInput if it's not this->GetInput().
  
    vtkUnstructuredGrid *redistributedInput = this->RedistributeDataSet(splitInput,
                                                                        input);
  
    if (redistributedInput == nullptr)
    {
      this->Kdtree->Delete();
      this->Kdtree = nullptr;
  
      vtkErrorMacro(<< "vtkDistributedDataFilter::Execute redistribute failure");
      return 1;
    }
  
    this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
  
    // Stage (3) - Add ghost cells to my sub grid.
  
    vtkUnstructuredGrid *expandedGrid = redistributedInput;
  
    if (this->GhostLevel > 0)
    {
      // Create global nodes IDs if we don't have them
  
      if (this->GetGlobalNodeIdArray(redistributedInput) == nullptr)
      {
        this->SetProgressText("Assign global point IDs");
        int rc = this->AssignGlobalNodeIds(redistributedInput);
        if (rc)
        {
          redistributedInput->Delete();
          this->Kdtree->Delete();
          this->Kdtree = nullptr;
          vtkErrorMacro(<< "vtkDistributedDataFilter::Execute global node id creation");
          return 1;
        }
      }
  
      // redistributedInput will be deleted by AcquireGhostCells
  
      this->SetProgressText("Exchange ghost cells");
      expandedGrid = this->AcquireGhostCells(redistributedInput);
    }
  
    // Stage (4) - Clip cells to the spatial region boundaries
  
    if (this->ClipCells)
    {
      this->SetProgressText("Clip boundary cells");
      this->ClipGridCells(expandedGrid);
      this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
    }
  
    // remove temporary arrays we created
  
    this->SetProgressText("Clean up and finish");
  
    vtkDataArray *da = expandedGrid->GetCellData()->GetArray(TEMP_ELEMENT_ID_NAME);
  
    if (da)
    {
      expandedGrid->GetCellData()->RemoveArray(TEMP_ELEMENT_ID_NAME);
    }
  
    da = expandedGrid->GetPointData()->GetArray(TEMP_NODE_ID_NAME);
  
    if (da)
    {
      expandedGrid->GetCellData()->RemoveArray(TEMP_NODE_ID_NAME);
    }
  
    output->ShallowCopy(expandedGrid);
    output->GetFieldData()->ShallowCopy(input->GetFieldData());
  
    expandedGrid->Delete();
  
    if (!this->RetainKdtree)
    {
      this->Kdtree->Delete();
      this->Kdtree = nullptr;
    }
    else
    {
      this->Kdtree->SetDataSet(nullptr);
    }
  
    this->UpdateProgress(1);
  
    return 1;
  }
  
  //----------------------------------------------------------------------------
  vtkUnstructuredGrid *vtkDistributedDataFilter::RedistributeDataSet(
    vtkDataSet *set, vtkDataSet *input)
  {
    // Create global cell ids before redistributing data.  These
    // will be necessary if we need ghost cells later on.
  
    vtkDataSet *inputPlus = set;
  
    if ((this->GhostLevel > 0) && (this->GetGlobalElementIdArray(set) == nullptr))
    {
      if (set == input)
      {
        inputPlus = set->NewInstance();
        inputPlus->ShallowCopy(set);
      }
  
      this->AssignGlobalElementIds(inputPlus);
    }
  
    // next call deletes inputPlus at the earliest opportunity
  
    vtkUnstructuredGrid *finalGrid = this->MPIRedistribute(inputPlus, input);
  
    return finalGrid;
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::PartitionDataAndAssignToProcesses(vtkDataSet *set)
  {
    if (this->Kdtree == nullptr)
    {
      this->Kdtree = vtkPKdTree::New();
      if (!this->UserCuts)
      {
        this->Kdtree->AssignRegionsContiguous();
      }
      this->Kdtree->SetTiming(this->GetTiming());
    }
    if (this->UserCuts)
    {
      this->Kdtree->SetCuts(this->UserCuts);
    }
  
    this->Kdtree->SetController(this->Controller);
    this->Kdtree->SetNumberOfRegionsOrMore(this->NumProcesses);
    this->Kdtree->SetMinCells(0);
    this->Kdtree->SetDataSet(set);
  
    // BuildLocator is smart enough to rebuild the k-d tree only if
    // the input geometry has changed, or the k-d tree build parameters
    // have changed.  It will reassign regions if the region assignment
    // scheme has changed.
  
    this->Kdtree->BuildLocator();
  
    int nregions = this->Kdtree->GetNumberOfRegions();
  
    if (nregions < this->NumProcesses)
    {
      if (nregions == 0)
      {
        vtkErrorMacro("Unable to build k-d tree structure");
      }
      else
      {
        vtkErrorMacro(<< "K-d tree must have at least one region per process.  "
                      << "Needed " << this->NumProcesses << ", has " << nregions);
      }
      this->Kdtree->Delete();
      this->Kdtree = nullptr;
      return 1;
    }
  
    if (this->Internals->UserRegionAssignments.size() > 0)
    {
      if (
        static_cast<int>(this->Internals->UserRegionAssignments.size()) !=
        nregions)
      {
        vtkWarningMacro("Mismatch in number of user-defined regions and regions"
          " the in KdTree. Ignoring user-defined regions.");
      }
      else
      {
        this->Kdtree->AssignRegions(
          &this->Internals->UserRegionAssignments[0], nregions);
      }
    }
  
    return 0;
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::ClipGridCells(vtkUnstructuredGrid *grid)
  {
    if (grid->GetNumberOfCells() == 0)
    {
      return 0;
    }
  
    // Global point IDs are meaningless after
    // clipping, since this tetrahedralizes the whole data set.
    // We remove that array.
  
    if (this->GetGlobalNodeIdArray(grid))
    {
      grid->GetPointData()->SetGlobalIds(nullptr);
    }
  
    this->ClipCellsToSpatialRegion(grid);
  
    return 0;
  }
  
  //----------------------------------------------------------------------------
  vtkUnstructuredGrid *
    vtkDistributedDataFilter::AcquireGhostCells(vtkUnstructuredGrid *grid)
  {
    if (this->GhostLevel < 1)
    {
      return grid;
    }
  
    // Create a search structure mapping global point IDs to local point IDs
  
    vtkIdType numPoints = grid->GetNumberOfPoints();
  
    vtkIdType *gnids = nullptr;
  
    if (numPoints > 0)
    {
      gnids = this->GetGlobalNodeIds(grid);
  
      if (!gnids)
      {
        vtkWarningMacro(<< "Can't create ghost cells without global node IDs");
        return grid;
      }
    }
  
    vtkDistributedDataFilterSTLCloak globalToLocalMap;
  
    for (int localPtId = 0; localPtId < numPoints; localPtId++)
    {
      const int id = gnids[localPtId];
      globalToLocalMap.IntMap.insert(std::pair<const int, int>(id, localPtId));
    }
  
    vtkUnstructuredGrid *expandedGrid= nullptr;
  
    if (this->IncludeAllIntersectingCells)
    {
      expandedGrid =
        this->AddGhostCellsDuplicateCellAssignment(grid, &globalToLocalMap);
    }
    else
    {
      expandedGrid =
        this->AddGhostCellsUniqueCellAssignment(grid, &globalToLocalMap);
    }
  
    convertGhostLevelsToBitFields(
      expandedGrid->GetCellData(), vtkDataSetAttributes::DUPLICATECELL);
    convertGhostLevelsToBitFields(
      expandedGrid->GetPointData(), vtkDataSetAttributes::DUPLICATEPOINT);
  
    return expandedGrid;
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::SingleProcessExecute(vtkDataSet *input,
                                                      vtkUnstructuredGrid *output)
  {
    vtkDebugMacro(<< "vtkDistributedDataFilter::SingleProcessExecute()");
  
    // we run the input through vtkMergeCells which will remove
    // duplicate points
  
    vtkDataSet* tmp = input->NewInstance();
    tmp->ShallowCopy(input);
  
    float tolerance = 0.0;
  
    if (this->RetainKdtree)
    {
      if (this->Kdtree == nullptr)
      {
        this->Kdtree = vtkPKdTree::New();
        if (this->UserCuts)
        {
          this->Kdtree->SetCuts(this->UserCuts);
        }
        this->Kdtree->SetTiming(this->GetTiming());
      }
  
      this->Kdtree->SetDataSet(tmp);
      this->Kdtree->BuildLocator();
      tolerance = (float)this->Kdtree->GetFudgeFactor();
      this->Kdtree->CreateGlobalDataArrayBounds();
    }
    else if (this->Kdtree)
    {
      this->Kdtree->Delete();
      this->Kdtree = nullptr;
    }
  
    vtkUnstructuredGrid *clean =
      vtkDistributedDataFilter::MergeGrids(&tmp, 1, DeleteYes,
                                           1, tolerance, 0);
  
    output->ShallowCopy(clean);
    clean->Delete();
  
    if (this->GhostLevel > 0)
    {
      // Add the vtkGhostType arrays.  We have the whole
      // data set, so all cells are level 0.
  
      vtkDistributedDataFilter::AddConstantUnsignedCharPointArray(
                                output, vtkDataSetAttributes::GhostArrayName(), 0);
      vtkDistributedDataFilter::AddConstantUnsignedCharCellArray(
                                output, vtkDataSetAttributes::GhostArrayName(), 0);
    }
  }
  
  //----------------------------------------------------------------------------
  void vtkDistributedDataFilter::ComputeMyRegionBounds()
  {
    delete [] this->ConvexSubRegionBounds;
    this->ConvexSubRegionBounds = nullptr;
  
    vtkIntArray *myRegions = vtkIntArray::New();
  
    this->Kdtree->GetRegionAssignmentList(this->MyId, myRegions);
  
    if (myRegions->GetNumberOfTuples() > 0)
    {
      this->NumConvexSubRegions =
        this->Kdtree->MinimalNumberOfConvexSubRegions(
          myRegions, &this->ConvexSubRegionBounds);
    }
    else
    {
      this->NumConvexSubRegions = 0;
    }
  
    myRegions->Delete();
  }
  
  //----------------------------------------------------------------------------
  int vtkDistributedDataFilter::CheckFieldArrayTypes(vtkDataSet *set)
  {
    int i;
  
    // problem - vtkIdType arrays are written out as int arrays
    // when marshalled with vtkDataWriter.  This is a problem
    // when receive the array and try to merge it with our own,
    // which is a vtkIdType
  
    vtkPointData *pd = set->GetPointData();
    vtkCellData *cd = set->GetCellData();
  
    int npointArrays = pd->GetNumberOfArrays();
  
    for (i=0; i<npointArrays; i++)
    {
      int arrayType = pd->GetArray(i)->GetDataType();
  
      if (arrayType == VTK_ID_TYPE)
      {
        return 1;
      }
    }
  
    int ncellArrays = cd->GetNumberOfArrays();
  
    for (i=0; i<ncellArrays; i++)
    {
      int arrayType = cd->GetArray(i)->GetDataType();
  
      if (arrayType == VTK_ID_TYPE)
      {
        return 1;
      }
    }
  
    return 0;
  }
  
  //-------------------------------------------------------------------------
  // Quickly spread input data around if there are more processes than
  // input data sets.
  //-------------------------------------------------------------------------
  struct vtkDistributedDataFilterProcInfo {
    vtkIdType had;
    int procId;
    vtkIdType has;
  };
  extern "C"
  {
    int vtkDistributedDataFilterSortSize(const void *s1, const void *s2)
    {
      vtkDistributedDataFilterProcInfo *a, *b;
  
      a = (struct vtkDistributedDataFilterProcInfo *)s1;
      b = (struct vtkDistributedDataFilterProcInfo *)s2;
  
      if (a->has < b->has)
      {
        return 1;
      }
      else if (a->has == b->has)
      {
        return 0;
      }
      else
      {
        return -1;
      }
    }
  }
  
  //----------------------------------------------------------------------------
  vtkDataSet *vtkDistributedDataFilter::TestFixTooFewInputFiles(vtkDataSet *input)
  {
    vtkIdType i;
    int proc;
    int me = this->MyId;
    int nprocs = this->NumProcesses;
  
    vtkIdType numMyCells = input->GetNumberOfCells();
  
    // Find out how many input cells each process has.
  
    vtkIdTypeArray *inputSize = this->ExchangeCounts(numMyCells, 0x0001);
    vtkIdType *sizes = inputSize->GetPointer(0);
  
    int *nodeType = new int [nprocs];
    const int Producer = 1;
    const int Consumer = 2;
    int numConsumers = 0;
    vtkIdType numTotalCells = 0;
  
    for (proc = 0; proc < nprocs ; proc++)
    {
      numTotalCells += sizes[proc];
      if (sizes[proc] == 0)
      {
        numConsumers++;
        nodeType[proc] = Consumer;
      }
      else
      {
        nodeType[proc] = Producer;
      }
    }
  
    if (numTotalCells == 0)
    {
      // Nothing to do.
      // Based on the comments in RequestData() where this method is called, if
      // this method returns nullptr, it indicates that there's no distribution to be
      // done. That's indeed the case for empty datasets. Hence we'll return nullptr.
      delete[] nodeType;
      inputSize->Delete();
      return nullptr;
    }
  
    if (numConsumers == 0)
    {
      // Nothing to do.  Every process has input data.
  
      delete [] nodeType;
      inputSize->Delete();
      return input;
    }
  
    vtkIdType cellsPerNode = numTotalCells / nprocs;
  
    vtkIdList **sendCells = new vtkIdList * [ nprocs ];
    memset(sendCells, 0, sizeof(vtkIdList *) * nprocs);
  
    if (numConsumers == nprocs - 1)
    {
      // Simple and common case.
      // Only one process has data and divides it among the rest.
  
      inputSize->Delete();
  
      if (nodeType[me] == Producer)
      {
        if (numTotalCells < nprocs)
        {
          // If there are not enough cells to go around, just give one cell
          // to each process, duplicating as necessary.
          for (proc = 0; proc < nprocs; proc++)
          {
            sendCells[proc] = vtkIdList::New();
            sendCells[proc]->SetNumberOfIds(1);
            sendCells[proc]->SetId(0, proc%numTotalCells);
          }
        }
        else
        {
          vtkIdType sizeLast = numTotalCells - ((nprocs-1) * cellsPerNode);
          vtkIdType cellId = 0;
  
          for (proc=0; proc<nprocs; proc++)
          {
            vtkIdType ncells = ((proc == nprocs - 1) ? sizeLast :cellsPerNode);
  
            sendCells[proc] = vtkIdList::New();
            sendCells[proc]->SetNumberOfIds(ncells);
  
            for (i=0; i<ncells; i++)
            {
              sendCells[proc]->SetId(i, cellId++);
            }
          }
        }
      }
    }
    else
    {
      if (numTotalCells < nprocs)
      {
        for (proc = 0; nodeType[proc] != Producer; proc++)
        {
          // empty loop.
        }
        if (proc == me)
        {
          // Have one process give out its cells to consumers.
          vtkIdType numCells = inputSize->GetValue(me);
          i = 0;
          sendCells[me] = vtkIdList::New();
          sendCells[me]->SetNumberOfIds(1);
          sendCells[me]->SetId(0, i++);
          if (i >= numCells) i = 0;
          for (proc = 0; proc < nprocs; proc++)
          {
            if (nodeType[proc] == Consumer)
            {
              sendCells[proc] = vtkIdList::New();
              sendCells[proc]->SetNumberOfIds(1);
              sendCells[proc]->SetId(0, i++);
              if (i >= numCells) i = 0;
            }
          }
        }
        else if (nodeType[me] == Producer)
        {
          // All other producers keep their own cells.
          vtkIdType numCells = inputSize->GetValue(me);
          sendCells[me] = vtkIdList::New();
          sendCells[me]->SetNumberOfIds(numCells);
          for (i = 0; i < numCells; i++)
          {
            sendCells[me]->SetId(i, i);
          }
        }
  
        inputSize->Delete();
      }
      else
      {
  
        // The processes with data send it to processes without data.
        // This is not the most balanced decomposition, and it is not the
        // fastest.  It is somewhere in between.
  
        vtkIdType minCells = (vtkIdType)(.8 * (double)cellsPerNode);
  
        struct vtkDistributedDataFilterProcInfo *procInfo
          = new struct vtkDistributedDataFilterProcInfo [nprocs];
  
        for (proc = 0; proc < nprocs ; proc++)
        {
          procInfo[proc].had   = inputSize->GetValue(proc);
          procInfo[proc].procId = proc;
          procInfo[proc].has   = inputSize->GetValue(proc);
        }
  
        inputSize->Delete();
  
        qsort(procInfo, nprocs, sizeof(struct vtkDistributedDataFilterProcInfo),
              vtkDistributedDataFilterSortSize);
  
        struct vtkDistributedDataFilterProcInfo *nextProducer = procInfo;
        struct vtkDistributedDataFilterProcInfo *nextConsumer
          = procInfo + (nprocs - 1);
  
        vtkIdType numTransferCells = 0;
  
        int sanityCheck=0;
        int nprocsSquared = nprocs * nprocs;
  
        while (sanityCheck++ < nprocsSquared)
        {
          int c = nextConsumer->procId;
  
          if (nodeType[c] == Producer)
          {
            break;
          }
  
          vtkIdType cGetMin = minCells - nextConsumer->has;
  
          if (cGetMin < 1)
          {
            nextConsumer--;
            continue;
          }
          vtkIdType cGetMax = cellsPerNode - nextConsumer->has;
  
          int p = nextProducer->procId;
  
          vtkIdType pSendMax = nextProducer->has - minCells;
  
          if (pSendMax < 1)
          {
            nextProducer++;
            continue;
          }
  
          vtkIdType transferSize = (pSendMax < cGetMax) ? pSendMax : cGetMax;
  
          if (me == p)
          {
            vtkIdType startCellId = nextProducer->had - nextProducer->has;
            sendCells[c] = vtkIdList::New();
            sendCells[c]->SetNumberOfIds(transferSize);
            for (i=0; i<transferSize; i++)
            {
              sendCells[c]->SetId(i, startCellId++);
            }
  
            numTransferCells += transferSize;
          }
  
          nextProducer->has -= transferSize;
          nextConsumer->has += transferSize;
  
          continue;
        }
  
        delete [] procInfo;
  
        if (sanityCheck > nprocsSquared)
        {
          vtkErrorMacro(<< "TestFixTooFewInputFiles error");
          for (i=0; i<nprocs; i++)
          {
            if (sendCells[i])
            {
              sendCells[i]->Delete();
            }
          }
          delete [] sendCells;
          delete [] nodeType;
          sendCells = nullptr;
        }
        else if (nodeType[me] == Producer)
        {
          vtkIdType keepCells = numMyCells - numTransferCells;
          vtkIdType startCellId = (vtkIdType)numTransferCells;
          sendCells[me] = vtkIdList::New();
          sendCells[me]->SetNumberOfIds(keepCells);
          for (i=0; i<keepCells; i++)
          {
            sendCells[me]->SetId(i, startCellId++);
          }
        }
      }
    }
  
    vtkUnstructuredGrid *newGrid = nullptr;
  
    if (sendCells)
    {
      newGrid = this->ExchangeMergeSubGrids(
             sendCells, DeleteYes, input, DeleteNo,
             DuplicateCellsNo, GhostCellsNo, 0x0011);
  
      delete [] sendCells;
      delete [] nodeType;
    }
  
    return newGrid;
  }
  
  //============================================================================
  // Communication routines - two versions:
  //   *Lean version use minimal memory
  //   *Fast versions use more memory, but are much faster
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::SetUpPairWiseExchange()
  {
    int iam = this->MyId;
    int nprocs = this->NumProcesses;
  
    delete [] this->Target;
    this->Target = nullptr;
  
    delete [] this->Source;
    this->Source = nullptr;
  
    if (nprocs == 1)
    {
      return;
    }
  
    this->Target = new int [nprocs - 1];
    this->Source = new int [nprocs - 1];
  
    for (int i=1; i< nprocs; i++)
    {
      this->Target[i-1] = (iam + i) % nprocs;
      this->Source[i-1] = (iam + nprocs - i) % nprocs;
    }
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::FreeIntArrays(vtkIdTypeArray **ar)
  {
    for (int i=0; i<this->NumProcesses; i++)
    {
      if (ar[i])
      {
        ar[i]->Delete();
      }
    }
  
    delete [] ar;
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::FreeIdLists(vtkIdList**lists, int nlists)
  {
    for (int i=0; i<nlists; i++)
    {
      if (lists[i])
      {
         lists[i]->Delete();
         lists[i] = nullptr;
      }
    }
  }
  
  //-------------------------------------------------------------------------
  vtkIdType vtkDistributedDataFilter::GetIdListSize(vtkIdList **lists, int nlists)
  {
    vtkIdType numCells = 0;
  
    for (int i=0; i<nlists; i++)
    {
      if (lists[i])
      {
        numCells += lists[i]->GetNumberOfIds();
      }
    }
  
    return numCells;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *
    vtkDistributedDataFilter::ExchangeMergeSubGrids(
                 vtkIdList **cellIds, int deleteCellIds,
                 vtkDataSet *myGrid, int deleteMyGrid,
                 int filterOutDuplicateCells, int ghostCellFlag,
                 int tag)
  {
    int nprocs = this->NumProcesses;
  
    int *numLists = new int [nprocs];
  
    vtkIdList ***listOfLists = new vtkIdList ** [nprocs];
  
    for (int i=0; i<nprocs; i++)
    {
      if (cellIds[i] == nullptr)
      {
        numLists[i] = 0;
      }
      else
      {
        numLists[i] = 1;
      }
  
      listOfLists[i] = &cellIds[i];
    }
  
    vtkUnstructuredGrid *grid = nullptr;
  
    if (this->UseMinimalMemory)
    {
      grid = this->ExchangeMergeSubGridsLean(listOfLists, numLists, deleteCellIds,
               myGrid, deleteMyGrid, filterOutDuplicateCells, ghostCellFlag, tag);
    }
    else
    {
      grid = this->ExchangeMergeSubGridsFast(listOfLists, numLists, deleteCellIds,
               myGrid, deleteMyGrid, filterOutDuplicateCells, ghostCellFlag, tag);
    }
  
    delete [] numLists;
    delete [] listOfLists;
  
    return grid;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *
    vtkDistributedDataFilter::ExchangeMergeSubGrids(
                 vtkIdList ***cellIds, int *numLists, int deleteCellIds,
                 vtkDataSet *myGrid, int deleteMyGrid,
                 int filterOutDuplicateCells, int ghostCellFlag,
                 int tag)
  {
    vtkUnstructuredGrid *grid = nullptr;
  
    if (this->UseMinimalMemory)
    {
      grid = this->ExchangeMergeSubGridsLean(cellIds, numLists, deleteCellIds,
               myGrid, deleteMyGrid, filterOutDuplicateCells, ghostCellFlag, tag);
    }
    else
    {
      grid = this->ExchangeMergeSubGridsFast(cellIds, numLists, deleteCellIds,
               myGrid, deleteMyGrid, filterOutDuplicateCells, ghostCellFlag, tag);
    }
    return grid;
  }
  
  //-------------------------------------------------------------------------
  vtkIdTypeArray *vtkDistributedDataFilter::ExchangeCounts(vtkIdType myCount, int tag)
  {
    vtkIdTypeArray *ia;
  
    if (this->UseMinimalMemory)
    {
      ia = this->ExchangeCountsLean(myCount, tag);
    }
    else
    {
      ia = this->ExchangeCountsFast(myCount, tag);
    }
    return ia;
  }
  
  //-------------------------------------------------------------------------
  vtkFloatArray **vtkDistributedDataFilter::
    ExchangeFloatArrays(vtkFloatArray **myArray, int deleteSendArrays, int tag)
  {
    vtkFloatArray **fa;
  
    if (this->UseMinimalMemory)
    {
      fa = this->ExchangeFloatArraysLean(myArray, deleteSendArrays, tag);
    }
    else
    {
      fa = this->ExchangeFloatArraysFast(myArray, deleteSendArrays, tag);
    }
    return fa;
  }
  
  //-------------------------------------------------------------------------
  vtkIdTypeArray **vtkDistributedDataFilter::
    ExchangeIdArrays(vtkIdTypeArray **myArray, int deleteSendArrays, int tag)
  {
    vtkIdTypeArray **ia;
  
    if (this->UseMinimalMemory)
    {
      ia = this->ExchangeIdArraysLean(myArray, deleteSendArrays, tag);
    }
    else
    {
      ia = this->ExchangeIdArraysFast(myArray, deleteSendArrays, tag);
    }
    return ia;
  }
  
  // ----------------------- Lean versions ----------------------------//
  vtkIdTypeArray *vtkDistributedDataFilter::ExchangeCountsLean(vtkIdType myCount, int tag)
  {
    vtkIdTypeArray *countArray = nullptr;
  
    vtkIdType i;
    int nprocs = this->NumProcesses;
  
    vtkMPICommunicator::Request req;
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    vtkIdType *counts = new vtkIdType [nprocs];
    counts[this->MyId] = myCount;
  
    if (!this->Source)
    {
      this->SetUpPairWiseExchange();
    }
  
    for (i = 0; i < this->NumProcesses - 1; i++)
    {
      int source = this->Source[i];
      int target = this->Target[i];
      mpiContr->NoBlockReceive(counts + source, 1, source, tag, req);
      mpiContr->Send(&myCount, 1, target, tag);
      req.Wait();
    }
  
    countArray = vtkIdTypeArray::New();
    countArray->SetArray(counts, nprocs, 0,
      vtkIdTypeArray::VTK_DATA_ARRAY_DELETE);
  
    return countArray;
  }
  
  //-------------------------------------------------------------------------
  vtkFloatArray **
    vtkDistributedDataFilter::ExchangeFloatArraysLean(vtkFloatArray **myArray,
                                                int deleteSendArrays, int tag)
  {
    vtkFloatArray **remoteArrays = nullptr;
  
    int i;
    int nprocs = this->NumProcesses;
    int me = this->MyId;
  
    vtkMPICommunicator::Request req;
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    int *recvSize = new int [nprocs];
    int *sendSize = new int [nprocs];
  
    if (!this->Source)
    {
      this->SetUpPairWiseExchange();
    }
  
    for (i= 0; i< nprocs; i++)
    {
      sendSize[i] = myArray[i] ? myArray[i]->GetNumberOfTuples() : 0;
      recvSize[i] = 0;
    }
  
    // Exchange sizes
  
    int nothers = nprocs - 1;
  
    for (i = 0; i < nothers; i++)
    {
      int source = this->Source[i];
      int target = this->Target[i];
      mpiContr->NoBlockReceive(recvSize + source, 1, source, tag, req);
      mpiContr->Send(sendSize + target, 1, target, tag);
      req.Wait();
    }
  
    // Exchange int arrays
  
    float **recvArrays = new float * [nprocs];
    memset(recvArrays, 0, sizeof(float *) * nprocs);
  
    if (sendSize[me] > 0)  // sent myself an array
    {
      recvSize[me] = sendSize[me];
      recvArrays[me] = new float [sendSize[me]];
      memcpy(recvArrays[me], myArray[me]->GetPointer(0), sendSize[me] * sizeof(float));
    }
  
    for (i = 0; i < nothers; i++)
    {
      int source = this->Source[i];
      int target = this->Target[i];
      recvArrays[source] = nullptr;
  
      if (recvSize[source] > 0)
      {
        recvArrays[source] = new float [recvSize[source]];
        if (recvArrays[source] == nullptr)
        {
          vtkErrorMacro(<<
            "vtkDistributedDataFilter::ExchangeIdArrays memory allocation");
          delete [] recvSize;
          delete [] recvArrays;
          return nullptr;
        }
        mpiContr->NoBlockReceive(recvArrays[source], recvSize[source], source, tag, req);
      }
  
      if (sendSize[target] > 0)
      {
        mpiContr->Send(myArray[target]->GetPointer(0), sendSize[target], target, tag);
      }
  
      if (myArray[target] && deleteSendArrays)
      {
        myArray[target]->Delete();
      }
  
      if (recvSize[source] > 0)
      {
        req.Wait();
      }
    }
  
    if (deleteSendArrays)
    {
      if (myArray[me])
      {
        myArray[me]->Delete();
      }
      delete [] myArray;
    }
  
    delete [] sendSize;
  
    remoteArrays = new vtkFloatArray * [nprocs];
  
    for (i=0; i<nprocs; i++)
    {
      if (recvSize[i] > 0)
      {
        remoteArrays[i] = vtkFloatArray::New();
        remoteArrays[i]->SetArray(recvArrays[i], recvSize[i],
          vtkFloatArray::VTK_DATA_ARRAY_DELETE);
      }
      else
      {
        remoteArrays[i] = nullptr;
      }
    }
  
    delete [] recvArrays;
    delete [] recvSize;
  
    return remoteArrays;
  }
  
  //-------------------------------------------------------------------------
  vtkIdTypeArray **
    vtkDistributedDataFilter::ExchangeIdArraysLean(vtkIdTypeArray **myArray,
                                                int deleteSendArrays, int tag)
  {
    vtkIdTypeArray **remoteArrays = nullptr;
  
    int i;
    int nprocs = this->NumProcesses;
    int me = this->MyId;
  
    vtkMPICommunicator::Request req;
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    vtkIdType *recvSize = new vtkIdType [nprocs];
    vtkIdType *sendSize = new vtkIdType [nprocs];
  
    if (!this->Source)
    {
      this->SetUpPairWiseExchange();
    }
  
    for (i= 0; i< nprocs; i++)
    {
      sendSize[i] = myArray[i] ? myArray[i]->GetNumberOfTuples() : 0;
      recvSize[i] = 0;
    }
  
    // Exchange sizes
  
    int nothers = nprocs - 1;
  
    for (i = 0; i < nothers; i++)
    {
      int source = this->Source[i];
      int target = this->Target[i];
      mpiContr->NoBlockReceive(recvSize + source, 1, source, tag, req);
      mpiContr->Send(sendSize + target, 1, target, tag);
      req.Wait();
    }
  
    // Exchange int arrays
  
    vtkIdType **recvArrays = new vtkIdType * [nprocs];
    memset(recvArrays, 0, sizeof(vtkIdType *) * nprocs);
  
    if (sendSize[me] > 0)  // sent myself an array
    {
      recvSize[me] = sendSize[me];
      recvArrays[me] = new vtkIdType [sendSize[me]];
      memcpy(recvArrays[me], myArray[me]->GetPointer(0), sendSize[me] * sizeof(vtkIdType));
    }
  
    for (i = 0; i < nothers; i++)
    {
      int source = this->Source[i];
      int target = this->Target[i];
      recvArrays[source] = nullptr;
  
      if (recvSize[source] > 0)
      {
        recvArrays[source] = new vtkIdType [recvSize[source]];
        if (recvArrays[source] == nullptr)
        {
          vtkErrorMacro(<<
            "vtkDistributedDataFilter::ExchangeIdArrays memory allocation");
          delete [] sendSize;
          delete [] recvSize;
          delete [] recvArrays;
          return nullptr;
        }
        mpiContr->NoBlockReceive(recvArrays[source], recvSize[source], source, tag, req);
      }
  
      if (sendSize[target] > 0)
      {
        mpiContr->Send(myArray[target]->GetPointer(0), sendSize[target], target, tag);
      }
  
      if (myArray[target] && deleteSendArrays)
      {
        myArray[target]->Delete();
      }
  
      if (recvSize[source] > 0)
      {
        req.Wait();
      }
    }
  
    if (deleteSendArrays)
    {
      if (myArray[me])
      {
        myArray[me]->Delete();
      }
      delete [] myArray;
    }
  
    delete [] sendSize;
  
    remoteArrays = new vtkIdTypeArray * [nprocs];
  
    for (i=0; i<nprocs; i++)
    {
      if (recvSize[i] > 0)
      {
        remoteArrays[i] = vtkIdTypeArray::New();
        remoteArrays[i]->SetArray(recvArrays[i], recvSize[i], 0,
          vtkIdTypeArray::VTK_DATA_ARRAY_DELETE);
      }
      else
      {
        remoteArrays[i] = nullptr;
      }
    }
  
    delete [] recvArrays;
    delete [] recvSize;
  
    return remoteArrays;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *
    vtkDistributedDataFilter::ExchangeMergeSubGridsLean(
      vtkIdList ***cellIds, int *numLists, int deleteCellIds,
      vtkDataSet *myGrid, int deleteMyGrid,
      int filterOutDuplicateCells,   // flag if different processes may send same cells
      int vtkNotUsed(ghostCellFlag),   // flag if these cells are ghost cells
      int tag)
  {
    vtkUnstructuredGrid *mergedGrid = nullptr;
    int i;
    int packedGridSendSize=0, packedGridRecvSize=0;
    char *packedGridSend=nullptr, *packedGridRecv=nullptr;
    int recvBufSize=0;
    int numReceivedGrids = 0;
  
    int nprocs = this->NumProcesses;
    int iam = this->MyId;
  
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
    vtkMPICommunicator::Request req;
  
    vtkDataSet *tmpGrid = myGrid->NewInstance();
    tmpGrid->ShallowCopy(myGrid);
  
    vtkDataSet **grids = new vtkDataSet * [nprocs];
  
    if (numLists[iam] > 0)
    {
      // I was extracting/packing/sending/unpacking ugrids of zero cells,
      // and this caused corrupted data structures.  I don't know why, but
      // I am now being careful not to do that.
  
      vtkIdType numCells =
        vtkDistributedDataFilter::GetIdListSize(cellIds[iam], numLists[iam]);
  
      if (numCells > 0)
      {
        grids[numReceivedGrids++] =
          this->ExtractCells(cellIds[iam], numLists[iam], deleteCellIds, tmpGrid);
      }
      else if (deleteCellIds)
      {
        vtkDistributedDataFilter::FreeIdLists(cellIds[iam], numLists[iam]);
      }
    }
  
    if (this->Source == nullptr)
    {
      this->SetUpPairWiseExchange();
    }
  
    int nothers = nprocs - 1;
  
    for (i=0; i<nothers; i++)
    {
      int target = this->Target[i];
      int source = this->Source[i];
  
      packedGridSendSize = 0;
  
      if (cellIds[target] && (numLists[target] > 0))
      {
        vtkIdType numCells = vtkDistributedDataFilter::GetIdListSize(
                        cellIds[target], numLists[target]);
  
        if (numCells > 0)
        {
          vtkUnstructuredGrid *sendGrid =
            this->ExtractCells(cellIds[target], numLists[target],
                                                 deleteCellIds, tmpGrid);
  
          packedGridSend = this->MarshallDataSet(sendGrid, packedGridSendSize);
          sendGrid->Delete();
        }
        else if (deleteCellIds)
        {
          vtkDistributedDataFilter::FreeIdLists(cellIds[target], numLists[target]);
        }
      }
  
      // exchange size of packed grids
  
      mpiContr->NoBlockReceive(&packedGridRecvSize, 1, source, tag, req);
      mpiContr->Send(&packedGridSendSize, 1, target, tag);
      req.Wait();
  
      if (packedGridRecvSize > recvBufSize)
      {
        delete [] packedGridRecv;
        packedGridRecv = new char [packedGridRecvSize];
        if (!packedGridRecv)
        {
          vtkErrorMacro(<<
            "vtkDistributedDataFilter::ExchangeMergeSubGrids memory allocation");
          delete [] grids;
          return nullptr;
        }
        recvBufSize = packedGridRecvSize;
      }
  
      if (packedGridRecvSize > 0)
      {
        mpiContr->NoBlockReceive(packedGridRecv, packedGridRecvSize, source, tag, req);
      }
  
      if (packedGridSendSize > 0)
      {
        mpiContr->Send(packedGridSend, packedGridSendSize, target, tag);
        delete [] packedGridSend;
      }
  
      if (packedGridRecvSize > 0)
      {
        req.Wait();
  
        grids[numReceivedGrids++] =
          this->UnMarshallDataSet(packedGridRecv, packedGridRecvSize);
      }
    }
  
    tmpGrid->Delete();
  
    if (recvBufSize > 0)
    {
      delete [] packedGridRecv;
      packedGridRecv = nullptr;
    }
  
    if (numReceivedGrids > 1)
    {
      // Merge received grids
  
      // this call will merge the grids and then delete them
  
      float tolerance = 0.0;
  
      if (this->Kdtree)
      {
        tolerance = (float)this->Kdtree->GetFudgeFactor();
      }
  
      mergedGrid =
        vtkDistributedDataFilter::MergeGrids(grids, numReceivedGrids, DeleteYes,
                                             1, tolerance, filterOutDuplicateCells);
  
    }
    else if (numReceivedGrids == 1)
    {
      mergedGrid = vtkUnstructuredGrid::SafeDownCast(grids[0]);
    }
    else
    {
      mergedGrid = this->ExtractZeroCellGrid(myGrid);
    }
  
    if (deleteMyGrid)
    {
      myGrid->Delete();
    }
  
    delete [] grids;
  
    return mergedGrid;
  }
  
  // ----------------------- Fast versions ----------------------------//
  vtkIdTypeArray *vtkDistributedDataFilter::ExchangeCountsFast(
    vtkIdType myCount, int vtkNotUsed(tag))
  {
    int nprocs = this->NumProcesses;
  
    vtkIdType *counts = new vtkIdType [nprocs];
    this->Controller->AllGather(&myCount, counts, 1);
  
    vtkIdTypeArray *countArray = vtkIdTypeArray::New();
    countArray->SetArray(counts, nprocs, 0,
      vtkIdTypeArray::VTK_DATA_ARRAY_DELETE);
    return countArray;
  }
  
  //-------------------------------------------------------------------------
  vtkFloatArray **
    vtkDistributedDataFilter::ExchangeFloatArraysFast(vtkFloatArray **myArray,
                                                int deleteSendArrays, int tag)
  {
    vtkFloatArray **fa = nullptr;
    int proc;
    int nprocs = this->NumProcesses;
    int iam = this->MyId;
  
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    int *sendSize = new int [nprocs];
    int *recvSize = new int [nprocs];
  
    for (proc=0; proc < nprocs; proc++)
    {
      recvSize[proc] = sendSize[proc] = 0;
  
      if (proc == iam)
      {
        continue;
      }
  
      if (myArray[proc])
      {
        sendSize[proc] = myArray[proc]->GetNumberOfTuples();
      }
    }
  
    // Exchange sizes of arrays to send and receive
  
    vtkMPICommunicator::Request *reqBuf = new vtkMPICommunicator::Request [nprocs];
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->NoBlockReceive(recvSize + proc, 1, proc, tag, reqBuf[proc]);
    }
  
    mpiContr->Barrier();
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->Send(sendSize + proc, 1, proc, tag);
    }
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      reqBuf[proc].Wait();
    }
  
    // Allocate buffers and post receives
  
    float **recvBufs = new float * [nprocs];
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (recvSize[proc] > 0)
      {
        recvBufs[proc] = new float [recvSize[proc]];
        mpiContr->NoBlockReceive(recvBufs[proc], recvSize[proc], proc, tag, reqBuf[proc]);
      }
      else
      {
        recvBufs[proc] = nullptr;
      }
    }
  
    mpiContr->Barrier();
  
    // Send all arrays
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (sendSize[proc] > 0)
      {
        mpiContr->Send(myArray[proc]->GetPointer(0), sendSize[proc], proc, tag);
      }
    }
    delete [] sendSize;
  
    // If I want to send an array to myself, place it in output now
  
    if (myArray[iam])
    {
      recvSize[iam] = myArray[iam]->GetNumberOfTuples();
      if (recvSize[iam] > 0)
      {
        recvBufs[iam] = new float [recvSize[iam]];
        memcpy(recvBufs[iam], myArray[iam]->GetPointer(0), recvSize[iam] * sizeof(float));
      }
    }
  
    if (deleteSendArrays)
    {
      for (proc=0; proc < nprocs; proc++)
      {
        if (myArray[proc])
        {
          myArray[proc]->Delete();
        }
      }
      delete [] myArray;
    }
  
    // Await incoming arrays
  
    fa = new vtkFloatArray * [nprocs];
    for (proc=0; proc < nprocs; proc++)
    {
      if (recvBufs[proc])
      {
        fa[proc] = vtkFloatArray::New();
        fa[proc]->SetArray(recvBufs[proc], recvSize[proc], 0,
          vtkFloatArray::VTK_DATA_ARRAY_DELETE);
      }
      else
      {
        fa[proc] = nullptr;
      }
    }
  
    delete [] recvSize;
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      if (recvBufs[proc])
      {
        reqBuf[proc].Wait();
      }
    }
  
    delete [] reqBuf;
    delete [] recvBufs;
  
    return fa;
  }
  
  //-------------------------------------------------------------------------
  vtkIdTypeArray **
    vtkDistributedDataFilter::ExchangeIdArraysFast(vtkIdTypeArray **myArray,
                                                int deleteSendArrays, int tag)
  {
    vtkIdTypeArray **ia = nullptr;
    int proc;
    int nprocs = this->NumProcesses;
    int iam = this->MyId;
  
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    vtkIdType *sendSize = new vtkIdType [nprocs];
    vtkIdType *recvSize = new vtkIdType [nprocs];
  
    for (proc=0; proc < nprocs; proc++)
    {
      recvSize[proc] = sendSize[proc] = 0;
  
      if (proc == iam)
      {
        continue;
      }
  
      if (myArray[proc])
      {
        sendSize[proc] = myArray[proc]->GetNumberOfTuples();
      }
    }
  
    // Exchange sizes of arrays to send and receive
  
    vtkMPICommunicator::Request *reqBuf = new vtkMPICommunicator::Request [nprocs];
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->NoBlockReceive(recvSize + proc, 1, proc, tag, reqBuf[proc]);
    }
  
    mpiContr->Barrier();
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->Send(sendSize + proc, 1, proc, tag);
    }
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      reqBuf[proc].Wait();
    }
  
    // Allocate buffers and post receives
  
    vtkIdType **recvBufs = new vtkIdType * [nprocs];
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (recvSize[proc] > 0)
      {
        recvBufs[proc] = new vtkIdType [recvSize[proc]];
        mpiContr->NoBlockReceive(recvBufs[proc], recvSize[proc], proc, tag, reqBuf[proc]);
      }
      else
      {
        recvBufs[proc] = nullptr;
      }
    }
  
    mpiContr->Barrier();
  
    // Send all arrays
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (sendSize[proc] > 0)
      {
        mpiContr->Send(myArray[proc]->GetPointer(0), sendSize[proc], proc, tag);
      }
    }
    delete [] sendSize;
  
    // If I want to send an array to myself, place it in output now
  
    if (myArray[iam])
    {
      recvSize[iam] = myArray[iam]->GetNumberOfTuples();
      if (recvSize[iam] > 0)
      {
        recvBufs[iam] = new vtkIdType [recvSize[iam]];
        memcpy(recvBufs[iam], myArray[iam]->GetPointer(0), recvSize[iam] * sizeof(vtkIdType));
      }
    }
  
    if (deleteSendArrays)
    {
      for (proc=0; proc < nprocs; proc++)
      {
        if (myArray[proc])
        {
          myArray[proc]->Delete();
        }
      }
      delete [] myArray;
    }
  
    // Await incoming arrays
  
    ia = new vtkIdTypeArray * [nprocs];
    for (proc=0; proc < nprocs; proc++)
    {
      if (recvBufs[proc])
      {
        ia[proc] = vtkIdTypeArray::New();
        ia[proc]->SetArray(recvBufs[proc], recvSize[proc], 0,
          vtkIdTypeArray::VTK_DATA_ARRAY_DELETE);
      }
      else
      {
        ia[proc] = nullptr;
      }
    }
  
    delete [] recvSize;
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      if (recvBufs[proc])
      {
        reqBuf[proc].Wait();
      }
    }
  
    delete [] reqBuf;
    delete [] recvBufs;
  
    return ia;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *
    vtkDistributedDataFilter::ExchangeMergeSubGridsFast(
      vtkIdList ***cellIds, int *numLists, int deleteCellIds,
      vtkDataSet *myGrid, int deleteMyGrid,
      int filterOutDuplicateCells,   // flag if different processes may send same cells
      int vtkNotUsed(ghostCellFlag),  // flag if these are ghost cells
      int tag)
  {
    vtkUnstructuredGrid *mergedGrid = nullptr;
    int proc;
    int nprocs = this->NumProcesses;
    int iam = this->MyId;
  
    vtkMPIController *mpiContr = vtkMPIController::SafeDownCast(this->Controller);
  
    vtkUnstructuredGrid **grids = new vtkUnstructuredGrid * [nprocs];
    char **sendBufs = new char * [nprocs];
    char **recvBufs = new char * [nprocs];
    int *sendSize = new int [nprocs];
    int *recvSize = new int [nprocs];
  
    // create & pack all sub grids
  
    vtkDataSet *tmpGrid = myGrid->NewInstance();
    tmpGrid->ShallowCopy(myGrid);
  
    for (proc=0; proc < nprocs; proc++)
    {
      recvSize[proc] = sendSize[proc] = 0;
      grids[proc] = nullptr;
      sendBufs[proc] = recvBufs[proc] = nullptr;
  
      if (numLists[proc] > 0)
      {
        vtkIdType numCells =
          vtkDistributedDataFilter::GetIdListSize(cellIds[proc], numLists[proc]);
  
        if (numCells > 0)
        {
          grids[proc] =
            vtkDistributedDataFilter::ExtractCells(cellIds[proc], numLists[proc],
                                            deleteCellIds, tmpGrid);
  
          if (proc != iam)
          {
            sendBufs[proc] = this->MarshallDataSet(grids[proc], sendSize[proc]);
            grids[proc]->Delete();
            grids[proc] = nullptr;
          }
        }
        else if (deleteCellIds)
        {
          vtkDistributedDataFilter::FreeIdLists(cellIds[proc], numLists[proc]);
        }
      }
    }
  
    tmpGrid->Delete();
  
    // Exchange sizes of grids to send and receive
  
    vtkMPICommunicator::Request *reqBuf = new vtkMPICommunicator::Request [nprocs];
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->NoBlockReceive(recvSize + proc, 1, proc, tag, reqBuf[proc]);
    }
  
    mpiContr->Barrier();
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      mpiContr->Send(sendSize + proc, 1, proc, tag);
    }
  
    for (proc=0; proc<nprocs; proc++)
    {
      if (proc == iam)
      {
        continue;
      }
      reqBuf[proc].Wait();
    }
  
    // Allocate buffers and post receives
  
    int numReceives = 0;
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (recvSize[proc] > 0)
      {
        recvBufs[proc] = new char [recvSize[proc]];
        mpiContr->NoBlockReceive(recvBufs[proc], recvSize[proc], proc, tag, reqBuf[proc]);
        numReceives++;
      }
    }
  
    mpiContr->Barrier();
  
    // Send all sub grids, then delete them
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (sendSize[proc] > 0)
      {
        mpiContr->Send(sendBufs[proc], sendSize[proc], proc, tag);
      }
    }
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (sendSize[proc] > 0)
      {
        delete [] sendBufs[proc];
      }
    }
  
    delete [] sendSize;
    delete [] sendBufs;
  
    // Await incoming sub grids, unpack them
  
    while (numReceives > 0)
    {
      for (proc=0; proc < nprocs; proc++)
      {
        if (recvBufs[proc] && (reqBuf[proc].Test() == 1))
        {
          grids[proc] = this->UnMarshallDataSet(recvBufs[proc], recvSize[proc]);
          delete [] recvBufs[proc];
          recvBufs[proc] = nullptr;
          numReceives--;
        }
      }
    }
  
    delete [] reqBuf;
    delete [] recvBufs;
    delete [] recvSize;
  
    // Merge received grids
  
    float tolerance = 0.0;
  
    if (this->Kdtree)
    {
      tolerance = (float)this->Kdtree->GetFudgeFactor();
    }
  
    int numReceivedGrids = 0;
  
    vtkDataSet **ds = new vtkDataSet * [nprocs];
  
    for (proc=0; proc < nprocs; proc++)
    {
      if (grids[proc] != nullptr)
      {
        ds[numReceivedGrids++] = static_cast<vtkDataSet *>(grids[proc]);
      }
    }
  
    delete [] grids;
  
    if (numReceivedGrids > 1)
    {
      // Normally, using this->GetGlobalNodeIds is the right thing.  However,
      // there is a bit of a bug here that this filter only works with ids
      // that are vtkIdType.  Otherwise, it will return nullptr as the global ids.
      // That is bad because then the global node ids will be stripped in the
      // MergeGrids method, and the number of point arrays will not match,
      // causing a crash latter on.
      //int useGlobalNodeIds = (this->GetGlobalNodeIds(ds[0]) != nullptr);
      int useGlobalNodeIds = (ds[0]->GetPointData()->GetGlobalIds() != nullptr);
  
      // this call will merge the grids and then delete them
      mergedGrid =
        vtkDistributedDataFilter::MergeGrids(ds, numReceivedGrids, DeleteYes,
                                             useGlobalNodeIds, tolerance,
                                             filterOutDuplicateCells);
  
    }
    else if (numReceivedGrids == 1)
    {
      mergedGrid = vtkUnstructuredGrid::SafeDownCast(ds[0]);
    }
    else
    {
      mergedGrid = this->ExtractZeroCellGrid(myGrid);
    }
  
    if (deleteMyGrid)
    {
      myGrid->Delete();
    }
  
    delete [] ds;
  
    return mergedGrid;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *vtkDistributedDataFilter::MPIRedistribute(vtkDataSet *in,
                                                                 vtkDataSet *input)
  {
    int proc;
    int nprocs = this->NumProcesses;
  
    // A cell belongs to a spatial region if it's centroid lies in that
    // region.  The kdtree object can create a list for each region of the
    // IDs of each cell I have read in that belong in that region.  If we
    // are building subgrids of all cells that intersect a region (a
    // superset of all cells that belong to a region) then the kdtree object
    // can build another set of lists of all cells that intersect each
    // region (but don't have their centroid in that region).
  
    if (this->IncludeAllIntersectingCells)
    {
      // TO DO:
      // We actually compute whether a cell intersects a spatial region.
      // This can be a lengthy calculation.  Perhaps it's good enough
      // to compute whether a cell's bounding box intersects the region.
      // Some of the cells we list will actually not be in the region, but
      // if we are clipping later, it doesn't matter.
      //
      // Is there any rendering algorithm that needs exactly all cells
      // which intersect the region, and no more?
  
      this->Kdtree->IncludeRegionBoundaryCellsOn();   // SLOW!!
    }
  
    this->Kdtree->CreateCellLists();  // required by GetCellIdsForProcess
  
    vtkIdList ***procCellLists = new vtkIdList ** [nprocs];
    int *numLists = new int [nprocs];
  
    for (proc = 0; proc < this->NumProcesses; proc++)
    {
      procCellLists[proc] = this->GetCellIdsForProcess(proc, numLists + proc);
    }
  
    int deleteDataSet = DeleteNo;
  
    if (in != input)
    {
      deleteDataSet = DeleteYes;
    }
  
    vtkUnstructuredGrid *myNewGrid =
      this->ExchangeMergeSubGrids(procCellLists, numLists, DeleteNo,
         in, deleteDataSet, DuplicateCellsNo, GhostCellsNo, 0x0012);
  
    for (proc = 0; proc < nprocs; proc++)
    {
      delete [] procCellLists[proc];
    }
  
    delete [] procCellLists;
    delete [] numLists;
  
    if (myNewGrid && (this->GhostLevel > 0))
    {
      vtkDistributedDataFilter::AddConstantUnsignedCharCellArray(
                              myNewGrid, vtkDataSetAttributes::GhostArrayName(), 0);
      vtkDistributedDataFilter::AddConstantUnsignedCharPointArray(
                              myNewGrid, vtkDataSetAttributes::GhostArrayName(), 0);
    }
    return myNewGrid;
  }
  
  //-------------------------------------------------------------------------
  char *vtkDistributedDataFilter::MarshallDataSet(vtkUnstructuredGrid *extractedGrid, int &len)
  {
    // taken from vtkCommunicator::WriteDataSet
  
    vtkUnstructuredGrid *copy;
    vtkDataSetWriter *writer = vtkDataSetWriter::New();
  
    copy = extractedGrid->NewInstance();
    copy->ShallowCopy(extractedGrid);
  
    // There is a problem with binary files with no data.
    if (copy->GetNumberOfCells() > 0)
    {
      writer->SetFileTypeToBinary();
    }
    writer->WriteToOutputStringOn();
    writer->SetInputData(copy);
  
    writer->Write();
  
    len = writer->GetOutputStringLength();
  
    char *packedFormat = writer->RegisterAndGetOutputString();
  
    writer->Delete();
  
    copy->Delete();
  
    return packedFormat;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid *vtkDistributedDataFilter::UnMarshallDataSet(char *buf, int size)
  {
    // taken from vtkCommunicator::ReadDataSet
  
    vtkDataSetReader *reader = vtkDataSetReader::New();
  
    reader->ReadFromInputStringOn();
  
    vtkCharArray* mystring = vtkCharArray::New();
  
    mystring->SetArray(buf, size, 1);
  
    reader->SetInputArray(mystring);
    mystring->Delete();
  
    vtkDataSet *output = reader->GetOutput();
    reader->Update();
  
    vtkUnstructuredGrid *newGrid = vtkUnstructuredGrid::New();
  
    newGrid->ShallowCopy(output);
  
    reader->Delete();
  
    return newGrid;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid
    *vtkDistributedDataFilter::ExtractCells(vtkIdList *cells, int deleteCellLists,
                                  vtkDataSet *in)
  {
    vtkIdList *tempCellList = nullptr;
  
    if (cells == nullptr)
    {
      // We'll get a zero cell unstructured grid which matches the input grid
      tempCellList = vtkIdList::New();
    }
    else
    {
      tempCellList = cells;
    }
  
    vtkUnstructuredGrid *subGrid = vtkDistributedDataFilter::ExtractCells(
      &tempCellList, 1, deleteCellLists, in);
  
    if (tempCellList != cells)
    {
      tempCellList->Delete();
    }
  
    return subGrid;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid
    *vtkDistributedDataFilter::ExtractCells(vtkIdList **cells, int nlists,
                   int deleteCellLists, vtkDataSet *in)
  {
    vtkDataSet* tmpInput = in->NewInstance();
    tmpInput->ShallowCopy(in);
  
    vtkExtractCells *extCells = vtkExtractCells::New();
  
    extCells->SetInputData(tmpInput);
  
    for (int i=0; i<nlists; i++)
    {
      if (cells[i])
      {
        extCells->AddCellList(cells[i]);
  
        if (deleteCellLists)
        {
          cells[i]->Delete();
        }
      }
    }
  
    extCells->Update();
  
    // If this process has no cells for these regions, a ugrid gets
    // created anyway with field array information
  
    vtkUnstructuredGrid *keepGrid = vtkUnstructuredGrid::New();
    keepGrid->ShallowCopy(extCells->GetOutput());
  
    extCells->Delete();
  
    tmpInput->Delete();
  
    return keepGrid;
  }
  
  //-------------------------------------------------------------------------
  vtkUnstructuredGrid
    *vtkDistributedDataFilter::ExtractZeroCellGrid(vtkDataSet *in)
  {
    vtkDataSet* tmpInput = in->NewInstance();
    tmpInput->ShallowCopy(in);
  
    vtkExtractCells *extCells = vtkExtractCells::New();
  
    extCells->SetInputData(tmpInput);
  
    extCells->Update();   // extract no cells
  
    vtkUnstructuredGrid *keepGrid = vtkUnstructuredGrid::New();
    keepGrid->ShallowCopy(extCells->GetOutput());
  
    extCells->Delete();
  
    tmpInput->Delete();
  
    return keepGrid;
  }
  
  //-------------------------------------------------------------------------
  
  // To save on storage, we return actual pointers into the vtkKdTree's lists
  // of cell IDs.  So don't free the memory they are pointing to.
  // vtkKdTree::DeleteCellLists will delete them all when we're done.
  
  vtkIdList **vtkDistributedDataFilter::GetCellIdsForProcess(int proc, int *nlists)
  {
    *nlists = 0;
  
    vtkIntArray *regions = vtkIntArray::New();
  
    int nregions = this->Kdtree->GetRegionAssignmentList(proc, regions);
  
    if (nregions == 0)
    {
      return nullptr;
    }
  
    *nlists = nregions;
  
    if (this->IncludeAllIntersectingCells)
    {
      *nlists *= 2;
    }
  
    vtkIdList **lists = new vtkIdList * [*nlists];
  
    int nextList = 0;
  
    for (int reg=0; reg < nregions; reg++)
    {
      lists[nextList++] = this->Kdtree->GetCellList(regions->GetValue(reg));
  
      if (this->IncludeAllIntersectingCells)
      {
        lists[nextList++] = this->Kdtree->GetBoundaryCellList(regions->GetValue(reg));
      }
    }
  
    regions->Delete();
  
    return lists;
  }
  
  //==========================================================================
  // Code related to clipping cells to the spatial region
  
  //-------------------------------------------------------------------------
  static int insideBoxFunction(vtkIdType cellId, vtkUnstructuredGrid *grid, void *data)
  {
    char *arrayName = (char *)data;
  
    vtkDataArray *da= grid->GetCellData()->GetArray(arrayName);
    vtkUnsignedCharArray *inside = vtkArrayDownCast<vtkUnsignedCharArray>(da);
  
    unsigned char where = inside->GetValue(cellId);
  
    return where;   // 1 if cell is inside spatial region, 0 otherwise
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::AddConstantUnsignedCharPointArray(
    vtkUnstructuredGrid *grid, const char *arrayName, unsigned char val)
  {
    vtkUnsignedCharArray *Array = vtkUnsignedCharArray::New();
    Array->SetName(arrayName);
  
    vtkIdType npoints = grid->GetNumberOfPoints();
    if (npoints)
    {
      unsigned char *vals = new unsigned char [npoints];
      memset(vals, val, npoints);
  
      Array->SetArray(vals, npoints, 0,
        vtkUnsignedCharArray::VTK_DATA_ARRAY_DELETE);
    }
  
    grid->GetPointData()->AddArray(Array);
  
    Array->Delete();
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::AddConstantUnsignedCharCellArray(
    vtkUnstructuredGrid *grid, const char *arrayName, unsigned char val)
  {
    vtkUnsignedCharArray *Array = vtkUnsignedCharArray::New();
    Array->SetName(arrayName);
  
    vtkIdType ncells = grid->GetNumberOfCells();
    if (ncells)
    {
      unsigned char *vals = new unsigned char [ncells];
      memset(vals, val, ncells);
  
      Array->SetArray(vals, ncells, 0, vtkUnsignedCharArray::VTK_DATA_ARRAY_DELETE);
    }
  
    grid->GetCellData()->AddArray(Array);
  
    Array->Delete();
  }
  
  #if 0
  //-------------------------------------------------------------------------
  // this is here temporarily, until vtkBoxClipDataSet is fixed to
  // be able to generate the clipped output
  void vtkDistributedDataFilter::ClipWithVtkClipDataSet(
             vtkUnstructuredGrid *grid, double *bounds,
             vtkUnstructuredGrid **outside, vtkUnstructuredGrid **inside)
  {
    vtkUnstructuredGrid *in;
    vtkUnstructuredGrid *out ;
  
    vtkClipDataSet *clipped = vtkClipDataSet::New();
  
    vtkBox *box = vtkBox::New();
    box->SetBounds(bounds);
  
    clipped->SetClipFunction(box);
    box->Delete();
    clipped->SetValue(0.0);
    clipped->InsideOutOn();
  
    clipped->SetInputData(grid);
  
    if (outside)
    {
      clipped->GenerateClippedOutputOn();
    }
  
    clipped->Update();
  
    if (outside)
    {
      out = clipped->GetClippedOutput();
      out->Register(this);
      *outside = out;
    }
  
    in = clipped->GetOutput();
    in->Register(this);
    *inside = in;
  
  
    clipped->Delete();
  }
  #endif
  
  //-------------------------------------------------------------------------
  // In general, vtkBoxClipDataSet is much faster and makes fewer errors.
  void vtkDistributedDataFilter::ClipWithBoxClipDataSet(
             vtkUnstructuredGrid *grid, double *bounds,
             vtkUnstructuredGrid **outside, vtkUnstructuredGrid **inside)
  {
    vtkUnstructuredGrid *in;
    vtkUnstructuredGrid *out ;
  
    vtkBoxClipDataSet *clipped = vtkBoxClipDataSet::New();
  
    clipped->SetBoxClip(bounds[0], bounds[1],
                        bounds[2], bounds[3], bounds[4], bounds[5]);
  
    clipped->SetInputData(grid);
  
    if (outside)
    {
      clipped->GenerateClippedOutputOn();
    }
  
    clipped->Update();
  
    if (outside)
    {
      out = clipped->GetClippedOutput();
      out->Register(this);
      *outside = out;
    }
  
    in = clipped->GetOutput();
    in->Register(this);
    *inside = in;
  
    clipped->Delete();
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::ClipCellsToSpatialRegion(vtkUnstructuredGrid *grid)
  {
    this->ComputeMyRegionBounds();
  
    if (this->NumConvexSubRegions > 1)
    {
      // here we would need to divide the grid into a separate grid for
      // each convex region, and then do the clipping
  
      vtkErrorMacro(<<
         "vtkDistributedDataFilter::ClipCellsToSpatialRegion - "
         "assigned regions do not form a single convex region");
  
      return ;
    }
  
    double *bounds = this->ConvexSubRegionBounds;
  
    if (this->GhostLevel > 0)
    {
      // We need cells outside the clip box as well.
  
      vtkUnstructuredGrid *outside;
      vtkUnstructuredGrid *inside;
  
  #if 1
      this->ClipWithBoxClipDataSet(grid, bounds, &outside, &inside);
  #else
      this->ClipWithVtkClipDataSet(grid, bounds, &outside, &inside);
  #endif
  
      grid->Initialize();
  
      // Mark the outside cells with a 0, the inside cells with a 1.
  
      int arrayNameLen = static_cast<int>(strlen(TEMP_INSIDE_BOX_FLAG));
      char *arrayName = new char [arrayNameLen + 1];
      strcpy(arrayName, TEMP_INSIDE_BOX_FLAG);
      vtkDistributedDataFilter::AddConstantUnsignedCharCellArray(outside, arrayName, 0);
      vtkDistributedDataFilter::AddConstantUnsignedCharCellArray(inside, arrayName, 1);
  
      // Combine inside and outside into a single ugrid.
  
      vtkDataSet *grids[2];
      grids[0] = inside;
      grids[1] = outside;
  
      vtkUnstructuredGrid *combined =
        vtkDistributedDataFilter::MergeGrids(grids, 2,  DeleteYes, 0,
                                             (float)this->Kdtree->GetFudgeFactor(), 0);
  
      // Extract the piece inside the box (level 0) and the requested
      // number of levels of ghost cells.
  
      vtkExtractUserDefinedPiece *ep = vtkExtractUserDefinedPiece::New();
  
      ep->SetConstantData(arrayName, arrayNameLen + 1);
      ep->SetPieceFunction(insideBoxFunction);
      ep->CreateGhostCellsOn();
  
      ep->GetExecutive()->GetOutputInformation(0)->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
        this->GhostLevel);
      ep->SetInputData(combined);
  
      ep->Update();
  
      grid->ShallowCopy(ep->GetOutput());
      grid->GetCellData()->RemoveArray(arrayName);
  
      ep->Delete();
      combined->Delete();
  
      delete [] arrayName;
    }
    else
    {
      vtkUnstructuredGrid *inside;
  
  #if 1
      this->ClipWithBoxClipDataSet(grid, bounds, nullptr, &inside);
  #else
      this->ClipWithVtkClipDataSet(grid, bounds, nullptr, &inside);
  #endif
  
      grid->ShallowCopy(inside);
      inside->Delete();
    }
  
    return;
  }
  
  //==========================================================================
  // Code related to assigning global node IDs and cell IDs
  
  //-------------------------------------------------------------------------
  int vtkDistributedDataFilter::AssignGlobalNodeIds(vtkUnstructuredGrid *grid)
  {
    int nprocs = this->NumProcesses;
    int pid;
    vtkIdType ptId;
    vtkIdType nGridPoints = grid->GetNumberOfPoints();
  
    vtkIdType *numPointsOutside = new vtkIdType [nprocs];
    memset(numPointsOutside, 0, sizeof(vtkIdType) * nprocs);
  
    vtkIdTypeArray *globalIds = vtkIdTypeArray::New();
    globalIds->SetNumberOfValues(nGridPoints);
    globalIds->SetName(TEMP_NODE_ID_NAME);
  
    // 1. Count the points in grid which lie within my assigned spatial region
  
    vtkIdType myNumPointsInside = 0;
  
    for (ptId = 0; ptId < nGridPoints; ptId++)
    {
      double *pt = grid->GetPoints()->GetPoint(ptId);
  
      if (this->InMySpatialRegion(pt[0], pt[1], pt[2]))
      {
        globalIds->SetValue(ptId, 0);  // flag it as mine
        myNumPointsInside++;
      }
      else
      {
        // Well, whose region is this point in?
  
        int regionId = this->Kdtree->GetRegionContainingPoint(pt[0],pt[1],pt[2]);
  
        pid = this->Kdtree->GetProcessAssignedToRegion(regionId);
  
        numPointsOutside[pid]++;
  
        pid += 1;
        pid *= -1;
  
        globalIds->SetValue(ptId, pid);  // a flag
      }
    }
  
    // 2. Gather and Broadcast this number of "Inside" points for each process.
  
    vtkIdTypeArray *numPointsInside = this->ExchangeCounts(myNumPointsInside, 0x0013);
  
    // 3. Assign global Ids to the points inside my spatial region
  
    vtkIdType firstId = 0;
    vtkIdType numGlobalIdsSoFar = 0;
  
    for (pid = 0; pid < nprocs; pid++)
    {
      if (pid < this->MyId)
      {
        firstId += numPointsInside->GetValue(pid);
      }
      numGlobalIdsSoFar += numPointsInside->GetValue(pid);
    }
  
    numPointsInside->Delete();
  
    for (ptId = 0; ptId < nGridPoints; ptId++)
    {
      if (globalIds->GetValue(ptId) == 0)
      {
        globalIds->SetValue(ptId, firstId++);
      }
    }
  
    // -----------------------------------------------------------------
    // All processes have assigned global IDs to the points in their grid
    // which lie within their assigned spatial region.
    // Now they have to get the IDs for the
    // points in their grid which lie outside their region, and which
    // are within the spatial region of another process.
    // -----------------------------------------------------------------
  
    // 4. For every other process, build a list of points I have
    // which are in the region of that process.  In practice, the
    // processes for which I need to request points IDs should be
    // a small subset of all the other processes.
  
    // question: if the vtkPointArray has type double, should we
    // send doubles instead of floats to insure we get the right
    // global ID back?
  
    vtkFloatArray **ptarrayOut = new vtkFloatArray * [nprocs];
    memset(ptarrayOut, 0, sizeof(vtkFloatArray *) * nprocs);
  
    vtkIdTypeArray **localIds     = new vtkIdTypeArray * [nprocs];
    memset(localIds, 0, sizeof(vtkIdTypeArray *) * nprocs);
  
    vtkIdType *next = new vtkIdType [nprocs];
    vtkIdType *next3 = new vtkIdType [nprocs];
  
    for (ptId = 0; ptId < nGridPoints; ptId++)
    {
      pid = globalIds->GetValue(ptId);
  
      if (pid >= 0)
      {
        continue;   // that's one of mine
      }
  
      pid *= -1;
      pid -= 1;
  
      if (ptarrayOut[pid] == nullptr)
      {
        vtkIdType npoints = numPointsOutside[pid];
  
        ptarrayOut[pid] = vtkFloatArray::New();
        ptarrayOut[pid]->SetNumberOfValues(npoints * 3);
  
        localIds[pid] = vtkIdTypeArray::New();
        localIds[pid]->SetNumberOfValues(npoints);
  
        next[pid] = 0;
        next3[pid] = 0;
      }
  
      localIds[pid]->SetValue(next[pid]++, ptId);
  
      double *dp = grid->GetPoints()->GetPoint(ptId);
  
      ptarrayOut[pid]->SetValue(next3[pid]++, (float)dp[0]);
      ptarrayOut[pid]->SetValue(next3[pid]++, (float)dp[1]);
      ptarrayOut[pid]->SetValue(next3[pid]++, (float)dp[2]);
    }
  
    delete [] numPointsOutside;
    delete [] next;
    delete [] next3;
  
    // 5. Do pairwise exchanges of the points we want global IDs for,
    //    and delete outgoing point arrays.
  
    vtkFloatArray **ptarrayIn = this->ExchangeFloatArrays(ptarrayOut,
                DeleteYes, 0x0014);
  
    // 6. Find the global point IDs that have been requested of me,
    //    and delete incoming point arrays.  Count "missing points":
    //    the number of unique points I receive which are not in my
    //    grid (this may happen if IncludeAllIntersectingCells is OFF).
  
    vtkIdType myNumMissingPoints = 0;
  
    vtkIdTypeArray **idarrayOut =
      this->FindGlobalPointIds(ptarrayIn, globalIds, grid, myNumMissingPoints);
  
    vtkIdTypeArray *missingCount = this->ExchangeCounts(myNumMissingPoints, 0x0015);
  
    if (this->IncludeAllIntersectingCells == 1)
    {
      // Make sure all points were found
  
      int aok = 1;
      for (pid=0; pid<nprocs; pid++)
      {
        if (missingCount->GetValue(pid) > 0)
        {
           vtkErrorMacro(<<
            "vtkDistributedDataFilter::AssignGlobalNodeIds bad point");
          aok = 0;
          break;
        }
      }
      if (!aok)
      {
        this->FreeIntArrays(idarrayOut);
        this->FreeIntArrays(localIds);
        missingCount->Delete();
        globalIds->Delete();
  
        return 1;
      }
    }
  
    // 7. Do pairwise exchanges of the global point IDs, and delete the
    //    outgoing point ID arrays.
  
    vtkIdTypeArray **idarrayIn = this->ExchangeIdArrays(idarrayOut,
                      DeleteYes, 0x0016);
  
    // 8. It's possible (if IncludeAllIntersectingCells is OFF) that some
    //    processes had "missing points".  Process A has a point P in it's
    //    grid which lies in the spatial region of process B.  But P is not
    //    in process B's grid.  We need to assign global IDs to these points
    //    too.
  
    vtkIdType *missingId = new vtkIdType [nprocs];
  
    if (this->IncludeAllIntersectingCells == 0)
    {
      missingId[0] = numGlobalIdsSoFar;
  
      for (pid = 1; pid < nprocs; pid++)
      {
        int prev = pid - 1;
        missingId[pid] = missingId[prev] + missingCount->GetValue(prev);
      }
    }
  
    missingCount->Delete();
  
    // 9. Update my ugrid with these mutually agreed upon global point IDs
  
    for (pid = 0; pid < nprocs; pid++)
    {
      if (idarrayIn[pid] == nullptr)
      {
        continue;
      }
  
      vtkIdType count = idarrayIn[pid]->GetNumberOfTuples();
  
      for (ptId = 0; ptId < count; ptId++)
      {
        vtkIdType myLocalId = localIds[pid]->GetValue(ptId);
        vtkIdType yourGlobalId = idarrayIn[pid]->GetValue(ptId);
  
        if (yourGlobalId >= 0)
        {
          globalIds->SetValue(myLocalId, yourGlobalId);
        }
        else
        {
          vtkIdType ptIdOffset = yourGlobalId * -1;
          ptIdOffset -= 1;
  
          globalIds->SetValue(myLocalId, missingId[pid] + ptIdOffset);
        }
      }
      localIds[pid]->Delete();
      idarrayIn[pid]->Delete();
    }
  
    delete [] localIds;
    delete [] idarrayIn;
    delete [] missingId;
  
    grid->GetPointData()->SetGlobalIds(globalIds);
    globalIds->Delete();
  
    return 0;
  }
  
  //-------------------------------------------------------------------------
  // If grids were distributed with IncludeAllIntersectingCells OFF, it's
  // possible there are points in my spatial region that are not in my
  // grid.  They need global Ids, so I will keep track of how many such unique
  // points I receive from other processes, and will assign them temporary
  // IDs.  They will get permanent IDs later on.
  
  vtkIdTypeArray **vtkDistributedDataFilter::FindGlobalPointIds(
       vtkFloatArray **ptarray, vtkIdTypeArray *ids, vtkUnstructuredGrid *grid,
       vtkIdType &numUniqueMissingPoints)
  {
    int nprocs = this->NumProcesses;
    vtkIdTypeArray **gids = new vtkIdTypeArray * [nprocs];
  
    if (grid->GetNumberOfCells() == 0)
    {
      // There are no cells in my assigned region
  
      memset(gids, 0, sizeof(vtkIdTypeArray *) * nprocs);
  
      return gids;
    }
  
    vtkKdTree *kd = vtkKdTree::New();
  
    kd->BuildLocatorFromPoints(grid->GetPoints());
  
    int procId;
    vtkIdType ptId, localId;
  
    vtkPointLocator *pl = nullptr;
    vtkPoints *missingPoints = nullptr;
  
    if (this->IncludeAllIntersectingCells == 0)
    {
      this->ComputeMyRegionBounds();
      pl = vtkPointLocator::New();
      pl->SetTolerance(this->Kdtree->GetFudgeFactor());
      missingPoints = vtkPoints::New();
      pl->InitPointInsertion(missingPoints, this->ConvexSubRegionBounds);
    }
  
    for (procId = 0; procId < nprocs; procId++)
    {
      if ((ptarray[procId] == nullptr) ||
          (ptarray[procId]->GetNumberOfTuples() == 0))
      {
        gids[procId] = nullptr;
        if (ptarray[procId]) ptarray[procId]->Delete();
        continue;
      }
  
      gids[procId] = vtkIdTypeArray::New();
  
      vtkIdType npoints = ptarray[procId]->GetNumberOfTuples() / 3;
  
      gids[procId]->SetNumberOfValues(npoints);
      vtkIdType next = 0;
  
      float *pt = ptarray[procId]->GetPointer(0);
  
      for (ptId = 0; ptId < npoints; ptId++)
      {
        localId = kd->FindPoint((double)pt[0], (double)pt[1], (double)pt[2]);
  
        if (localId >= 0)
        {
          gids[procId]->SetValue(next++, ids->GetValue(localId));  // global Id
        }
        else
        {
          // This point is not in my grid
  
          if (this->IncludeAllIntersectingCells)
          {
            // This is an error
            gids[procId]->SetValue(next++, -1);
            numUniqueMissingPoints++;
          }
          else
          {
            // Flag these with a negative point ID.  We'll assign
            // them real point IDs later.
  
            vtkIdType nextId;
            double dpt[3];
            dpt[0] = pt[0]; dpt[1] = pt[1]; dpt[2] = pt[2];
            pl->InsertUniquePoint(dpt, nextId);
  
            nextId += 1;
            nextId *= -1;
            gids[procId]->SetValue(next++, nextId);
          }
        }
        pt += 3;
      }
  
      ptarray[procId]->Delete();
    }
  
    delete [] ptarray;
  
    kd->Delete();
  
    if (missingPoints)
    {
      numUniqueMissingPoints = missingPoints->GetNumberOfPoints();
      missingPoints->Delete();
      pl->Delete();
    }
  
    return gids;
  }
  
  //-------------------------------------------------------------------------
  int vtkDistributedDataFilter::AssignGlobalElementIds(vtkDataSet *in)
  {
    vtkIdType i;
    vtkIdType myNumCells = in->GetNumberOfCells();
    vtkIdTypeArray *numCells = this->ExchangeCounts(myNumCells, 0x0017);
  
    vtkIdTypeArray *globalCellIds = vtkIdTypeArray::New();
    globalCellIds->SetNumberOfValues(myNumCells);
    //DDM - do we need to mark this as the GID array?
    globalCellIds->SetName(TEMP_ELEMENT_ID_NAME);
  
    vtkIdType StartId = 0;
  
    for (i=0; i < this->MyId; i++)
    {
      StartId += numCells->GetValue(i);
    }
  
    numCells->Delete();
  
    for (i=0; i<myNumCells; i++)
    {
      globalCellIds->SetValue(i, StartId++);
    }
  
    in->GetCellData()->SetGlobalIds(globalCellIds);
  
    globalCellIds->Delete();
  
    return 0;
  }
  
  //========================================================================
  // Code related to acquiring ghost cells
  
  //-------------------------------------------------------------------------
  int vtkDistributedDataFilter::InMySpatialRegion(float x, float y, float z)
  {
    return this->InMySpatialRegion((double)x, (double)y, (double)z);
  }
  int vtkDistributedDataFilter::InMySpatialRegion(double x, double y, double z)
  {
    this->ComputeMyRegionBounds();
  
    double *box = this->ConvexSubRegionBounds;
  
    if (!box)
    {
      return 0;
    }
  
    // To avoid ambiguity, a point on a boundary is assigned to
    // the region for which it is on the upper boundary.  Or
    // (in one dimension) the region between points A and B
    // contains all points p such that A < p <= B.
  
    if ( (x <= box[0]) || (x > box[1]) ||
         (y <= box[2]) || (y > box[3]) ||
         (z <= box[4]) || (z > box[5])   )
    {
        return 0;
    }
  
    return 1;
  }
  
  //-----------------------------------------------------------------------
  int vtkDistributedDataFilter::StrictlyInsideMyBounds(float x, float y, float z)
  {
    return this->StrictlyInsideMyBounds((double)x, (double)y, (double)z);
  }
  
  //-----------------------------------------------------------------------
  int vtkDistributedDataFilter::StrictlyInsideMyBounds(double x, double y, double z)
  {
    this->ComputeMyRegionBounds();
  
    double *box = this->ConvexSubRegionBounds;
  
    if (!box)
    {
      return 0;
    }
  
    if ( (x <= box[0]) || (x >= box[1]) ||
         (y <= box[2]) || (y >= box[3]) ||
         (z <= box[4]) || (z >= box[5])   )
    {
        return 0;
    }
  
    return 1;
  }
  
  //-----------------------------------------------------------------------
  vtkIdTypeArray **vtkDistributedDataFilter::MakeProcessLists(
                                      vtkIdTypeArray **pointIds,
                                      vtkDistributedDataFilterSTLCloak *procs)
  {
    // Build a list of pointId/processId pairs for each process that
    // sent me point IDs.  The process Ids are all those processes
    // that had the specified point in their ghost level zero grid.
  
    int nprocs = this->NumProcesses;
  
    std::multimap<int, int>::iterator mapIt;
  
    vtkIdTypeArray **processList = new vtkIdTypeArray * [nprocs];
    memset(processList, 0, sizeof (vtkIdTypeArray *) * nprocs);
  
    for (int i=0; i<nprocs; i++)
    {
      if (pointIds[i] == nullptr)
      {
        continue;
      }
  
      vtkIdType size = pointIds[i]->GetNumberOfTuples();
  
      if (size > 0)
      {
        for (vtkIdType j=0; j<size; )
        {
          // These are all the points in my spatial region
          // for which process "i" needs ghost cells.
  
          vtkIdType gid = pointIds[i]->GetValue(j);
          vtkIdType ncells = pointIds[i]->GetValue(j+1);
  
          mapIt = procs->IntMultiMap.find(gid);
  
          while (mapIt != procs->IntMultiMap.end() && mapIt->first == gid)
          {
            int processId = mapIt->second;
  
            if (processId != i)
            {
              // Process "i" needs to know that process
              // "processId" also has cells using this point
  
              if (processList[i] == nullptr)
              {
                processList[i] = vtkIdTypeArray::New();
              }
              processList[i]->InsertNextValue(gid);
              processList[i]->InsertNextValue(processId);
            }
            ++mapIt;
          }
          j += (2 + ncells);
        }
      }
    }
  
    return processList;
  }
  
  //-----------------------------------------------------------------------
  vtkIdTypeArray *vtkDistributedDataFilter::AddPointAndCells(
                          vtkIdType gid, vtkIdType localId, vtkUnstructuredGrid *grid,
                          vtkIdType *gidCells, vtkIdTypeArray *ids)
  {
    if (ids == nullptr)
    {
      ids = vtkIdTypeArray::New();
    }
  
    ids->InsertNextValue(gid);
  
    vtkIdList *cellList = vtkIdList::New();
  
    grid->GetPointCells(localId, cellList);
  
    vtkIdType numCells = cellList->GetNumberOfIds();
  
    ids->InsertNextValue(numCells);
  
    for (vtkIdType j=0; j<numCells; j++)
    {
      vtkIdType globalCellId = gidCells[cellList->GetId(j)];
      ids->InsertNextValue(globalCellId);
    }
  
    cellList->Delete();
  
    return ids;
  }
  
  //-----------------------------------------------------------------------
  vtkIdTypeArray **vtkDistributedDataFilter::GetGhostPointIds(
    int ghostLevel, vtkUnstructuredGrid *grid,
    int AddCellsIAlreadyHave)
  {
    int nprocs = this->NumProcesses;
    int me = this->MyId;
    vtkIdType numPoints = grid->GetNumberOfPoints();
  
    vtkIdTypeArray **ghostPtIds = new vtkIdTypeArray * [nprocs];
    memset(ghostPtIds, 0, sizeof(vtkIdTypeArray *) * nprocs);
  
    if (numPoints < 1)
    {
      return ghostPtIds;
    }
  
    int processId = -1;
    int regionId = -1;
  
    vtkPKdTree *kd = this->Kdtree;
  
    vtkPoints *pts = grid->GetPoints();
  
    vtkIdType *gidsPoint = this->GetGlobalNodeIds(grid);
    vtkIdType *gidsCell = this->GetGlobalElementIds(grid);
  
  
    vtkUnsignedCharArray *uca = grid->GetPointGhostArray();
    unsigned char *levels = uca->GetPointer(0);
  
    unsigned char level = (unsigned char)(ghostLevel - 1);
  
    for (vtkIdType i=0; i<numPoints; i++)
    {
      double *pt = pts->GetPoint(i);
      regionId = kd->GetRegionContainingPoint(pt[0], pt[1], pt[2]);
      processId = kd->GetProcessAssignedToRegion(regionId);
  
      if (ghostLevel == 1)
      {
        // I want all points that are outside my spatial region
  
        if (processId == me)
        {
          continue;
        }
  
        // Don't include points that are not part of any cell
  
        int used = vtkDistributedDataFilter::LocalPointIdIsUsed(grid, i);
  
        if (!used)
        {
          continue;
        }
      }
      else if (levels[i] != level)
      {
        continue; // I want all points having the correct ghost level
      }
  
      vtkIdType gid = gidsPoint[i];
  
      if (AddCellsIAlreadyHave)
      {
        // To speed up exchange of ghost cells and creation of
        // new ghost cell grid, we tell other
        // processes which cells we already have, so they don't
        // send them to us.
  
        ghostPtIds[processId] =
          vtkDistributedDataFilter::AddPointAndCells(gid, i, grid, gidsCell,
                                         ghostPtIds[processId]);
      }
      else
      {
        if (ghostPtIds[processId] == nullptr)
        {
          ghostPtIds[processId] = vtkIdTypeArray::New();
        }
        ghostPtIds[processId]->InsertNextValue(gid);
        ghostPtIds[processId]->InsertNextValue(0);
      }
    }
    return ghostPtIds;
  }
  
  //-----------------------------------------------------------------------
  int vtkDistributedDataFilter::LocalPointIdIsUsed(
                                vtkUnstructuredGrid *grid, int ptId)
  {
    int used = 1;
  
    int numPoints = grid->GetNumberOfPoints();
  
    if ((ptId < 0) || (ptId >= numPoints))
    {
      used = 0;
    }
    else
    {
      vtkIdType id = (vtkIdType)ptId;
      vtkIdList *cellList = vtkIdList::New();
  
      grid->GetPointCells(id, cellList);
  
      if (cellList->GetNumberOfIds() == 0)
      {
        used = 0;
      }
  
      cellList->Delete();
    }
  
    return used;
  }
  
  //-----------------------------------------------------------------------
  int vtkDistributedDataFilter::GlobalPointIdIsUsed(vtkUnstructuredGrid *grid,
                      int ptId, vtkDistributedDataFilterSTLCloak *globalToLocal)
  {
    int used = 1;
  
    std::map<int, int>::iterator mapIt;
  
    mapIt = globalToLocal->IntMap.find(ptId);
  
    if (mapIt == globalToLocal->IntMap.end())
    {
      used = 0;
    }
    else
    {
      int id = mapIt->second;
  
      used = vtkDistributedDataFilter::LocalPointIdIsUsed(grid, id);
    }
  
    return used;
  }
  
  //-----------------------------------------------------------------------
  vtkIdType vtkDistributedDataFilter::FindId(vtkIdTypeArray *ids, vtkIdType gid, vtkIdType startLoc)
  {
    vtkIdType gidLoc = -1;
  
    if (ids == nullptr)
    {
      return gidLoc;
    }
  
    vtkIdType numIds = ids->GetNumberOfTuples();
  
    while ((ids->GetValue(startLoc) != gid) && (startLoc < numIds))
    {
      vtkIdType ncells = ids->GetValue(++startLoc);
      startLoc += (ncells + 1);
    }
  
    if (startLoc < numIds)
    {
      gidLoc = startLoc;
    }
  
    return gidLoc;
  }
  
  //-----------------------------------------------------------------------
  // We create an expanded grid with the required number of ghost
  // cells.  This is for the case where IncludeAllIntersectingCells is OFF.
  // This means that when the grid was redistributed, each cell was
  // uniquely assigned to one process, the process owning the spatial
  // region that the cell's centroid lies in.
  vtkUnstructuredGrid *
  vtkDistributedDataFilter::AddGhostCellsUniqueCellAssignment(
                                 vtkUnstructuredGrid *myGrid,
                                 vtkDistributedDataFilterSTLCloak *globalToLocalMap)
  {
    int i,j,k;
    int ncells=0;
    int processId=0;
    int gid=0;
    vtkIdType size=0;
  
    int nprocs = this->NumProcesses;
    int me = this->MyId;
  
    int gl = 1;
  
    // For each ghost level, processes request and send ghost cells
  
    vtkUnstructuredGrid *newGhostCellGrid = nullptr;
    vtkIdTypeArray **ghostPointIds = nullptr;
  
    vtkDistributedDataFilterSTLCloak *insidePointMap =
      new vtkDistributedDataFilterSTLCloak;
    std::multimap<int, int>::iterator mapIt;
  
    while (gl <= this->GhostLevel)
    {
      // For ghost level 1, create a list for each process (not
      // including me) of all points I have in that process'
      // assigned region.  We use this list for two purposes:
      // (1) to build a list on each process of all other processes
      // that have cells containing points in our region, (2)
      // these are some of the points that we need ghost cells for.
      //
      // For ghost level above 1, create a list for each process
      // (including me) of all my points in that process' assigned
      // region for which I need ghost cells.
  
      if (gl == 1)
      {
        ghostPointIds = this->GetGhostPointIds(gl, myGrid, 0);
      }
      else
      {
        ghostPointIds = this->GetGhostPointIds(gl, newGhostCellGrid, 1);
      }
  
      // Exchange these lists.
  
      vtkIdTypeArray **insideIds =
        this->ExchangeIdArrays(ghostPointIds, DeleteNo,
                                0x0018);
  
      if (gl == 1)
      {
        // For every point in my region that was sent to me by another process,
        // I now know the identity of all processes having cells containing
        // that point.  Begin by building a mapping from point IDs to the IDs
        // of processes that sent me that point.
  
        for (i=0; i<nprocs; i++)
        {
          if (insideIds[i] == nullptr)
          {
            continue;
          }
  
          size = insideIds[i]->GetNumberOfTuples();
  
          if (size > 0)
          {
            for (j=0; j<size; j+=2)
            {
              // map global point id to process ids
              const int id = (int)insideIds[i]->GetValue(j);
              insidePointMap->IntMultiMap.insert(std::pair<const int, int>(id, i));
            }
          }
        }
      }
  
      // Build a list of pointId/processId pairs for each process that
      // sent me point IDs.  To process P, for every point ID sent to me
      // by P, I send the ID of every other process (not including myself
      // and P) that has cells in it's ghost level 0 grid which use
      // this point.
  
      vtkIdTypeArray **processListSent
        = this->MakeProcessLists(insideIds, insidePointMap);
  
      // Exchange these new lists.
  
      vtkIdTypeArray **processList =
        this->ExchangeIdArrays(processListSent, DeleteYes,
                                0x0019);
  
      // I now know the identity of every process having cells containing
      // points I need ghost cells for.  Create a request to each process
      // for these cells.
  
      vtkIdTypeArray **ghostCellsPlease = new vtkIdTypeArray * [nprocs];
      for (i=0; i<nprocs; i++)
      {
        ghostCellsPlease[i] = vtkIdTypeArray::New();
        ghostCellsPlease[i]->SetNumberOfComponents(1);
      }
  
      for (i=0; i<nprocs; i++)
      {
        if (i == me)
        {
          continue;
        }
  
        if (ghostPointIds[i])       // points I have in your spatial region,
        {                         // maybe you have cells that use them?
  
          for (j=0; j<ghostPointIds[i]->GetNumberOfTuples(); j++)
          {
            ghostCellsPlease[i]->InsertNextValue(ghostPointIds[i]->GetValue(j));
          }
        }
        if (processList[i])         // other processes you say that also have
        {                         // cells using those points
          size = processList[i]->GetNumberOfTuples();
          vtkIdType *array = processList[i]->GetPointer(0);
          int nextLoc = 0;
  
          for (j=0; j < size; j += 2)
          {
            gid = array[j];
            processId = array[j+1];
  
            ghostCellsPlease[processId]->InsertNextValue(gid);
  
            if (gl > 1)
            {
              // add the list of cells I already have for this point
  
              int where =
                vtkDistributedDataFilter::FindId(ghostPointIds[i], gid, nextLoc);
  
              if (where < 0)
              {
                // error really, not sure what to do
                nextLoc = 0;
                ghostCellsPlease[processId]->InsertNextValue(0);
                continue;
              }
  
              ncells = ghostPointIds[i]->GetValue(where + 1);
  
              ghostCellsPlease[processId]->InsertNextValue(ncells);
  
              for (k=0; k <ncells; k++)
              {
                vtkIdType cellId = ghostPointIds[i]->GetValue(where + 2 + k);
                ghostCellsPlease[processId]->InsertNextValue(cellId);
              }
  
              nextLoc = where;
            }
            else
            {
              ghostCellsPlease[processId]->InsertNextValue(0);
            }
          }
        }
        if ((gl==1) && insideIds[i])   // points you have in my spatial region,
        {                            // which I may need ghost cells for
          for (j=0; j<insideIds[i]->GetNumberOfTuples();)
          {
            gid = insideIds[i]->GetValue(j);
            int used = vtkDistributedDataFilter::GlobalPointIdIsUsed(
                                    myGrid, gid, globalToLocalMap);
            if (used)
            {
              ghostCellsPlease[i]->InsertNextValue(gid);
              ghostCellsPlease[i]->InsertNextValue(0);
            }
  
            ncells = insideIds[i]->GetValue(j+1);
            j += (ncells + 2);
          }
        }
      }
  
      if (gl > 1)
      {
        if (ghostPointIds[me])   // these points are actually inside my region
        {
          size = ghostPointIds[me]->GetNumberOfTuples();
  
          for (i=0; i<size;)
          {
            gid = ghostPointIds[me]->GetValue(i);
            ncells = ghostPointIds[me]->GetValue(i+1);
  
            mapIt = insidePointMap->IntMultiMap.find(gid);
  
            if (mapIt != insidePointMap->IntMultiMap.end())
            {
              while (mapIt->first == gid)
              {
                processId = mapIt->second;
                ghostCellsPlease[processId]->InsertNextValue(gid);
                ghostCellsPlease[processId]->InsertNextValue(ncells);
  
                for (k=0; k<ncells; k++)
                {
                  vtkIdType cellId = ghostPointIds[me]->GetValue(i+1+k);
                  ghostCellsPlease[processId]->InsertNextValue(cellId);
                }
  
                ++mapIt;
              }
            }
            i += (ncells + 2);
          }
        }
      }
  
      this->FreeIntArrays(ghostPointIds);
      this->FreeIntArrays(insideIds);
      this->FreeIntArrays(processList);
  
      // Exchange these ghost cell requests.
  
      vtkIdTypeArray **ghostCellRequest
        = this->ExchangeIdArrays(ghostCellsPlease, DeleteYes,
                                  0x001a);
  
      // Build a list of cell IDs satisfying each request received.
      // Delete request arrays.
  
      vtkIdList **sendCellList =
        this->BuildRequestedGrids(ghostCellRequest, myGrid, globalToLocalMap);
  
      // Build subgrids and exchange them
  
      vtkUnstructuredGrid *incomingGhostCells = this->ExchangeMergeSubGrids(
               sendCellList, DeleteYes, myGrid, DeleteNo, DuplicateCellsNo,
               GhostCellsYes, 0x001b);
  
      delete [] sendCellList;
  
      // Set ghost level of new cells, and merge into grid of other
      // ghost cells received.
  
      newGhostCellGrid = this->SetMergeGhostGrid(newGhostCellGrid,
                                incomingGhostCells, gl, globalToLocalMap);
  
      this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
  
      gl++;
    }
  
    delete insidePointMap;
  
    vtkUnstructuredGrid *newGrid = nullptr;
  
    if (newGhostCellGrid && (newGhostCellGrid->GetNumberOfCells() > 0))
    {
      vtkDataSet *grids[2];
  
      grids[0] = myGrid;
      grids[1] = newGhostCellGrid;
  
      int useGlobalNodeIds = (this->GetGlobalNodeIds(myGrid)?1:0);
  
      newGrid =
        vtkDistributedDataFilter::MergeGrids(grids, 2, DeleteYes, useGlobalNodeIds, 0, 0);
    }
    else
    {
      newGrid = myGrid;
    }
  
    return newGrid;
  }
  
  //-----------------------------------------------------------------------
  // We create an expanded grid that contains the ghost cells we need.
  // This is in the case where IncludeAllIntersectingCells is ON.  This
  // is easier in some respects because we know if that if a point lies
  // in a region owned by a particular process, that process has all
  // cells which use that point.  So it is easy to find ghost cells.
  // On the otherhand, because cells are not uniquely assigned to regions,
  // we may get multiple processes sending us the same cell, so we
  // need to filter these out.
  vtkUnstructuredGrid *
  vtkDistributedDataFilter::AddGhostCellsDuplicateCellAssignment(
                               vtkUnstructuredGrid *myGrid,
                               vtkDistributedDataFilterSTLCloak *globalToLocalMap)
  {
    int i,j;
  
    int nprocs = this->NumProcesses;
    int me = this->MyId;
  
    int gl = 1;
  
    // For each ghost level, processes request and send ghost cells
  
    vtkUnstructuredGrid *newGhostCellGrid = nullptr;
    vtkIdTypeArray **ghostPointIds = nullptr;
    vtkIdTypeArray **extraGhostPointIds = nullptr;
  
    std::map<int, int>::iterator mapIt;
  
    vtkPoints *pts = myGrid->GetPoints();
  
    while (gl <= this->GhostLevel)
    {
      // For ghost level 1, create a list for each process of points
      // in my grid which lie in that other process' spatial region.
      // This is normally all the points for which I need ghost cells,
      // with one EXCEPTION.  If a cell is axis-aligned, and a face of
      // the cell is on my upper boundary, then the vertices of this
      // face are in my spatial region, but I need their ghost cells.
      // I can detect this case when the process across the boundary
      // sends me a request for ghost cells of these points.
      //
      // For ghost level above 1, create a list for each process of
      // points in my ghost grid which are in that process' spatial
      // region and for which I need ghost cells.
  
      if (gl == 1)
      {
        ghostPointIds = this->GetGhostPointIds(gl, myGrid, 1);
      }
      else
      {
        ghostPointIds = this->GetGhostPointIds(gl, newGhostCellGrid, 1);
      }
  
      // Exchange these lists.
  
      vtkIdTypeArray **insideIds =
        this->ExchangeIdArrays(ghostPointIds, DeleteYes, 0x001c);
  
      // For ghost level 1, examine the points Ids I received from
      // other processes, to see if the exception described above
      // applies and I need ghost cells from them for those points.
  
      if (gl == 1)
      {
        vtkIdType *gidsCell = this->GetGlobalElementIds(myGrid);
  
        extraGhostPointIds = new vtkIdTypeArray * [nprocs];
  
        for (i=0; i<nprocs; i++)
        {
          extraGhostPointIds[i] = nullptr;
  
          if (i == me)
          {
            continue;
          }
  
          if (insideIds[i] == nullptr)
          {
            continue;
          }
  
          vtkIdType size = insideIds[i]->GetNumberOfTuples();
  
          for (j=0; j<size;)
          {
            vtkIdType gid = insideIds[i]->GetValue(j);
            vtkIdType ncells = insideIds[i]->GetValue(j+1);
            j += (ncells + 2);
  
            mapIt = globalToLocalMap->IntMap.find(gid);
  
            if (mapIt == globalToLocalMap->IntMap.end())
            {
              // This point must be right on my boundary, and
              // not connected to any cell intersecting my region.
  
              continue;
            }
            vtkIdType localId = mapIt->second;
  
            double *pt = pts->GetPoint(localId);
  
            int interior = this->StrictlyInsideMyBounds(pt[0], pt[1], pt[2]);
  
            if (!interior)
            {
              extraGhostPointIds[i] = this->AddPointAndCells(gid, localId,
                              myGrid, gidsCell, extraGhostPointIds[i]);
            }
          }
        }
  
        // Exchange these lists.
  
        vtkIdTypeArray **extraInsideIds =
          this->ExchangeIdArrays(extraGhostPointIds, DeleteYes, 0x001d);
  
        // Add the extra point ids to the previous list
  
        for (i=0; i<nprocs; i++)
        {
          if (i == me)
          {
            continue;
          }
  
          if (extraInsideIds[i])
          {
            vtkIdType size = extraInsideIds[i]->GetNumberOfTuples();
  
            if (insideIds[i] == nullptr)
            {
              insideIds[i] = vtkIdTypeArray::New();
            }
  
            for (j=0; j<size; j++)
            {
              insideIds[i]->InsertNextValue(extraInsideIds[i]->GetValue(j));
            }
          }
        }
          this->FreeIntArrays(extraInsideIds);
      }
  
      // Build a list of cell IDs satisfying each request received.
  
      vtkIdList **sendCellList =
        this->BuildRequestedGrids(insideIds, myGrid, globalToLocalMap);
  
      // Build subgrids and exchange them
  
      vtkUnstructuredGrid *incomingGhostCells =
        this->ExchangeMergeSubGrids( sendCellList, DeleteYes, myGrid, DeleteNo,
                                     DuplicateCellsYes, GhostCellsYes, 0x001e);
  
      delete [] sendCellList;
  
      // Set ghost level of new cells, and merge into grid of other
      // ghost cells received.
  
      newGhostCellGrid = this->SetMergeGhostGrid(newGhostCellGrid,
                                incomingGhostCells, gl, globalToLocalMap);
  
      this->UpdateProgress(this->NextProgressStep++ * this->ProgressIncrement);
  
      gl++;
    }
  
    vtkUnstructuredGrid *newGrid = nullptr;
  
    if (newGhostCellGrid && (newGhostCellGrid->GetNumberOfCells() > 0))
    {
      vtkDataSet *grids[2];
  
      grids[0] = myGrid;
      grids[1] = newGhostCellGrid;
  
      int useGlobalNodeIds = (this->GetGlobalNodeIds(myGrid)?1:0);
      newGrid =
        vtkDistributedDataFilter::MergeGrids(grids, 2, DeleteYes,
                                             useGlobalNodeIds, 0, 0);
    }
    else
    {
      newGrid = myGrid;
    }
  
    return newGrid;
  }
  
  //-----------------------------------------------------------------------
  // For every process that sent me a list of point IDs, create a list
  // of all the cells I have in my original grid containing those points.
  // We omit cells the remote process already has.
  
  vtkIdList **vtkDistributedDataFilter::BuildRequestedGrids(
                          vtkIdTypeArray **globalPtIds,
                          vtkUnstructuredGrid *grid,
                          vtkDistributedDataFilterSTLCloak *ptIdMap)
  {
    vtkIdType id;
    int proc;
    int nprocs = this->NumProcesses;
    vtkIdType cellId;
    vtkIdType nelts;
  
    // for each process, create a list of the ids of cells I need
    // to send to it
  
    std::map<int, int>::iterator imap;
  
    vtkIdList *cellList = vtkIdList::New();
  
    vtkIdList **sendCells = new vtkIdList * [nprocs];
  
    for (proc = 0; proc < nprocs; proc++)
    {
      sendCells[proc] = vtkIdList::New();
  
      if (globalPtIds[proc] == nullptr)
      {
        continue;
      }
  
      if ((nelts = globalPtIds[proc]->GetNumberOfTuples()) == 0)
      {
        globalPtIds[proc]->Delete();
        continue;
      }
  
      vtkIdType *ptarray = globalPtIds[proc]->GetPointer(0);
  
      std::set<vtkIdType> subGridCellIds;
  
      vtkIdType nYourCells = 0;
  
      for (id = 0; id < nelts; id += (nYourCells + 2))
      {
        vtkIdType ptId = ptarray[id];
        nYourCells = ptarray[id+1];
  
        imap = ptIdMap->IntMap.find(ptId);
  
        if (imap == ptIdMap->IntMap.end())
        {
          continue; // I don't have this point
        }
  
        vtkIdType myPtId = (vtkIdType)imap->second;   // convert to my local point Id
  
        grid->GetPointCells(myPtId, cellList);
  
        vtkIdType nMyCells = cellList->GetNumberOfIds();
  
        if (nMyCells == 0)
        {
          continue;
        }
  
        if (nYourCells > 0)
        {
          // We don't send cells the remote process tells us it already
          // has.  This is much faster than removing duplicate cells on
          // the receive side.
  
          vtkIdType *remoteCells = ptarray + id + 2;
          vtkIdType *gidCells = this->GetGlobalElementIds(grid);
  
          vtkDistributedDataFilter::RemoveRemoteCellsFromList(cellList,
                                       gidCells, remoteCells, nYourCells);
        }
  
        vtkIdType nSendCells = cellList->GetNumberOfIds();
  
        if (nSendCells == 0)
        {
          continue;
        }
  
        for (cellId = 0; cellId < nSendCells; cellId++)
        {
          subGridCellIds.insert(cellList->GetId(cellId));
        }
      }
  
      globalPtIds[proc]->Delete();
  
      vtkIdType numUniqueCellIds = static_cast<vtkIdType>(subGridCellIds.size());
  
      if (numUniqueCellIds == 0)
      {
        continue;
      }
  
      sendCells[proc]->SetNumberOfIds(numUniqueCellIds);
      vtkIdType next = 0;
  
      std::set<vtkIdType>::iterator it;
  
      for (it = subGridCellIds.begin(); it != subGridCellIds.end(); ++it)
      {
        sendCells[proc]->SetId(next++, *it);
      }
    }
  
    delete [] globalPtIds;
  
    cellList->Delete();
  
    return sendCells;
  }
  
  //-----------------------------------------------------------------------
  void vtkDistributedDataFilter::RemoveRemoteCellsFromList(
    vtkIdList *cellList, vtkIdType *gidCells, vtkIdType *remoteCells, vtkIdType nRemoteCells)
  {
    vtkIdType id, nextId;
    vtkIdType id2;
    vtkIdType nLocalCells = cellList->GetNumberOfIds();
  
    // both lists should be very small, so we just do an n^2 lookup
  
    for (id = 0, nextId = 0; id < nLocalCells; id++)
    {
      vtkIdType localCellId  = cellList->GetId(id);
      vtkIdType globalCellId = gidCells[localCellId];
  
      int found = 0;
  
      for (id2 = 0; id2 < nRemoteCells; id2++)
      {
        if (remoteCells[id2] == globalCellId)
        {
          found = 1;
          break;
        }
      }
  
      if (!found)
      {
        cellList->SetId(nextId++, localCellId);
      }
    }
  
    cellList->SetNumberOfIds(nextId);
  }
  
  //-----------------------------------------------------------------------
  // Set the ghost levels for the points and cells in the received cells.
  // Merge the new ghost cells into the supplied grid, and return the new grid.
  // Delete all grids except the new merged grid.
  
  vtkUnstructuredGrid *vtkDistributedDataFilter::SetMergeGhostGrid(
                              vtkUnstructuredGrid *ghostCellGrid,
                              vtkUnstructuredGrid *incomingGhostCells,
                              int ghostLevel, vtkDistributedDataFilterSTLCloak *idMap)
  
  {
    int i;
  
    if (incomingGhostCells->GetNumberOfCells() < 1)
    {
      return ghostCellGrid;
    }
  
    // Set the ghost level of all new cells, and set the ghost level of all
    // the points.  We know some points in the new grids actually have ghost
    // level one lower, because they were on the boundary of the previous
    // grid.  This is OK if ghostLevel is > 1.  When we merge, vtkMergeCells
    // will skip these points because they are already in the previous grid.
    // But if ghostLevel is 1, those boundary points were in our original
    // grid, and we need to use the global ID map to determine if the
    // point ghost levels should be set to 0.
  
    vtkUnsignedCharArray *cellGL = incomingGhostCells->GetCellGhostArray();
    vtkUnsignedCharArray *ptGL  = incomingGhostCells->GetPointGhostArray();
  
  
    unsigned char *ia = cellGL->GetPointer(0);
  
    for (i=0; i < incomingGhostCells->GetNumberOfCells(); i++)
    {
      ia[i] = (unsigned char)ghostLevel;
    }
  
    ia = ptGL->GetPointer(0);
  
    for (i=0; i < incomingGhostCells->GetNumberOfPoints(); i++)
    {
      ia[i] = (unsigned char)ghostLevel;
    }
  
    // now merge
  
    vtkUnstructuredGrid *mergedGrid = incomingGhostCells;
  
    if (ghostCellGrid && (ghostCellGrid->GetNumberOfCells() > 0))
    {
      vtkDataSet *sets[2];
  
      sets[0] = ghostCellGrid;     // both sets will be deleted by MergeGrids
      sets[1] = incomingGhostCells;
  
      int useGlobalNodeIds = (this->GetGlobalNodeIds(ghostCellGrid)?1:0);
      mergedGrid =
        vtkDistributedDataFilter::MergeGrids(sets, 2, DeleteYes, useGlobalNodeIds, 0.0, 0);
    }
  
    // If this is ghost level 1, mark any points from our original grid
    // as ghost level 0.
  
    if (ghostLevel == 1)
    {
      ptGL = mergedGrid->GetPointGhostArray();
  
      vtkIdType *gidPoints = this->GetGlobalNodeIds(mergedGrid);
      int npoints = mergedGrid->GetNumberOfPoints();
  
      std::map<int, int>::iterator imap;
  
      for (i=0; i < npoints; i++)
      {
        imap = idMap->IntMap.find(gidPoints[i]);
  
        if (imap != idMap->IntMap.end())
        {
          ptGL->SetValue(i,0);   // found among my ghost level 0 cells
        }
      }
    }
  
    return mergedGrid;
  }
  
  //-----------------------------------------------------------------------
  vtkUnstructuredGrid *vtkDistributedDataFilter::MergeGrids(
    vtkDataSet **sets, int nsets, int deleteDataSets,
    int useGlobalNodeIds, float pointMergeTolerance,
    int useGlobalCellIds)
  {
    int i;
  
    if (nsets == 0)
    {
      return nullptr;
    }
  
    vtkUnstructuredGrid *newGrid = vtkUnstructuredGrid::New();
    // Any global ids should be consistent, so make sure they are passed.
    newGrid->GetPointData()->CopyGlobalIdsOn();
    newGrid->GetCellData()->CopyGlobalIdsOn();
  
    vtkMergeCells *mc = vtkMergeCells::New();
    mc->SetUnstructuredGrid(newGrid);
  
    mc->SetTotalNumberOfDataSets(nsets);
  
    vtkIdType totalPoints = 0;
    vtkIdType totalCells = 0;
  
    for (i=0; i<nsets; i++)
    {
      totalPoints += sets[i]->GetNumberOfPoints();
      totalCells += sets[i]->GetNumberOfCells();
      // Only use global ids if they are available.
      useGlobalNodeIds = (   useGlobalNodeIds
                          && (sets[i]->GetPointData()->GetGlobalIds() != nullptr) );
      useGlobalCellIds = (   useGlobalCellIds
                          && (sets[i]->GetCellData()->GetGlobalIds() != nullptr) );
    }
  
    mc->SetTotalNumberOfPoints(totalPoints);
    mc->SetTotalNumberOfCells(totalCells);
  
    if (!useGlobalNodeIds)
    {
      mc->SetPointMergeTolerance(pointMergeTolerance);
    }
    mc->SetUseGlobalIds(useGlobalNodeIds);
    mc->SetUseGlobalCellIds(useGlobalCellIds);
  
    for (i=0; i<nsets; i++)
    {
      mc->MergeDataSet(sets[i]);
  
      if (deleteDataSets)
      {
        sets[i]->Delete();
      }
    }
  
    mc->Finish();
    mc->Delete();
  
    return newGrid;
  }
  
  //-------------------------------------------------------------------------
  int vtkDistributedDataFilter::RequestDataObject(vtkInformation*,
    vtkInformationVector** inputVector,
    vtkInformationVector* outputVector)
  {
    vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
    if (!inInfo)
    {
      return 0;
    }
  
    vtkDataObject *input = vtkDataObject::GetData(inInfo);
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    if (input)
    {
      vtkDataObject *output = vtkDataObject::GetData(outInfo);
      // If input is composite dataset, output is a vtkMultiBlockDataSet of
      // unstructrued grids.
      // If input is a dataset, output is an unstructured grid.
      if (!output ||
        (input->IsA("vtkCompositeDataSet") && !output->IsA("vtkMultiBlockDataSet")) ||
        (input->IsA("vtkDataSet") && !output->IsA("vtkUnstructuredGrid")))
      {
        vtkDataObject* newOutput = nullptr;
        if (input->IsA("vtkCompositeDataSet"))
        {
          newOutput = vtkMultiBlockDataSet::New();
        }
        else // if (input->IsA("vtkDataSet"))
        {
          newOutput = vtkUnstructuredGrid::New();
        }
        outInfo->Set(vtkDataObject::DATA_OBJECT(), newOutput);
        newOutput->Delete();
      }
      return 1;
    }
  
    return 0;
  }
  
  //-------------------------------------------------------------------------
  int vtkDistributedDataFilter::FillInputPortInformation(int, vtkInformation *info)
  {
    info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkCompositeDataSet");
    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
    return 1;
  }
  
  //-------------------------------------------------------------------------
  vtkPKdTree *vtkDistributedDataFilter::GetKdtree()
  {
    if (this->Kdtree == nullptr)
    {
      this->Kdtree = vtkPKdTree::New();
      this->Kdtree->AssignRegionsContiguous();
      this->Kdtree->SetTiming(this->GetTiming());
    }
  
    return this->Kdtree;
  }
  
  //-------------------------------------------------------------------------
  void vtkDistributedDataFilter::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
  
    os << indent << "Kdtree: " << this->Kdtree << endl;
    os << indent << "Controller: " << this->Controller << endl;
    os << indent << "NumProcesses: " << this->NumProcesses << endl;
    os << indent << "MyId: " << this->MyId << endl;
    os << indent << "Target: " << this->Target << endl;
    os << indent << "Source: " << this->Source << endl;
    os << indent << "RetainKdtree: " << this->RetainKdtree << endl;
    os << indent << "IncludeAllIntersectingCells: " << this->IncludeAllIntersectingCells << endl;
    os << indent << "ClipCells: " << this->ClipCells << endl;
  
    os << indent << "Timing: " << this->Timing << endl;
    os << indent << "UseMinimalMemory: " << this->UseMinimalMemory << endl;
  }
  