// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen // SPDX-License-Identifier: BSD-3-Clause #include "vtkGenerateGlobalIds.h" #include "vtkBoundingBox.h" #include "vtkCellData.h" #include "vtkCompositeDataIterator.h" #include "vtkCompositeDataSet.h" #include "vtkDIYExplicitAssigner.h" #include "vtkDIYUtilities.h" #include "vtkDataSet.h" #include "vtkIdList.h" #include "vtkIdTypeArray.h" #include "vtkLogger.h" #include "vtkMath.h" #include "vtkMultiProcessController.h" #include "vtkNew.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" #include "vtkPoints.h" #include "vtkSMPThreadLocalObject.h" #include "vtkSMPTools.h" #include "vtkStaticPointLocator.h" #include "vtkTuple.h" #include "vtkUnsignedCharArray.h" #include "vtkUnstructuredGrid.h" #include #include #include #include // clang-format off #include "vtk_diy2.h" // #define DIY_USE_SPDLOG #include VTK_DIY2(diy/mpi.hpp) #include VTK_DIY2(diy/master.hpp) #include VTK_DIY2(diy/link.hpp) #include VTK_DIY2(diy/reduce.hpp) #include VTK_DIY2(diy/reduce-operations.hpp) #include VTK_DIY2(diy/partners/swap.hpp) #include VTK_DIY2(diy/assigner.hpp) #include VTK_DIY2(diy/algorithms.hpp) // clang-format on namespace impl { VTK_ABI_NAMESPACE_BEGIN static vtkBoundingBox AllReduceBounds( diy::mpi::communicator& comm, std::vector> points) { vtkBoundingBox bbox; for (auto& pts : points) { if (pts) { double bds[6]; pts->GetBounds(bds); bbox.AddBounds(bds); } } vtkDIYUtilities::AllReduce(comm, bbox); return bbox; } /** * This is the main implementation of the global id generation algorithm. * The code is similar for both point and cell ids generation except small * differences that are implemented using `BlockTraits` and `BlockT`. * * The general algorithm can be described as: * - sort points (or cells) globally so that all "coincident" points (or cells) * are within the same block; * - merge coincident points (or cells) per block and assign unique ids for * unique points (or cells) -- note this is local to each block since we know * all coincident points are same block after earlier step; * - uniquify the generated ids globally by exchanging information of local * unique id counts; * - communicate back the assigned unique id to the source block where the point * (or cell) came from. */ template static bool GenerateIds(vtkDataObject* dobj, vtkGenerateGlobalIds* self, bool cell_centers) { const double tolerance = self->GetTolerance(); self->UpdateProgress(0.0); diy::mpi::communicator comm = vtkDIYUtilities::GetCommunicator(self->GetController()); vtkLogStartScope(TRACE, "extract points"); auto datasets = vtkCompositeDataSet::GetDataSets(dobj); datasets.erase(std::remove_if(datasets.begin(), datasets.end(), [&cell_centers](vtkDataSet* ds) { return ds == nullptr || ds->GetNumberOfPoints() == 0 || (cell_centers && ds->GetNumberOfCells() == 0); }), datasets.end()); const auto points = vtkDIYUtilities::ExtractPoints(datasets, cell_centers); vtkLogEndScope("extract points"); // get the bounds for the domain globally. const diy::ContinuousBounds gdomain = vtkDIYUtilities::Convert(impl::AllReduceBounds(comm, points)); const int local_num_blocks = static_cast(points.size()); vtkDIYExplicitAssigner assigner(comm, local_num_blocks, /*pow-of-2*/ true); diy::Master master( comm, 1, -1, []() { return static_cast(new ElementBlockT); }, [](void* b) -> void { delete static_cast(b); }); vtkLogStartScope(TRACE, "populate master"); std::vector gids; assigner.local_gids(comm.rank(), gids); for (size_t lid = 0; lid < gids.size(); ++lid) { auto block = new ElementBlockT(); if (lid < points.size() && points[lid] != nullptr) { assert(datasets[lid] != nullptr); block->Initialize(gids[lid], points[lid], datasets[lid], cell_centers); } auto link = new diy::RegularContinuousLink(3, gdomain, gdomain); master.add(gids[lid], block, link); } vtkLogEndScope("populate master"); self->UpdateProgress(0.25); if (assigner.nblocks() > 1) { vtkLogStartScope(TRACE, "kdtree"); // use diy::kdtree to shuffle points around so that all spatially co-located // points are within a block. diy::kdtree(master, assigner, 3, gdomain, &ElementBlockT::Elements, /*hist_bins=*/512); vtkLogEndScope("kdtree"); } self->UpdateProgress(0.50); vtkLogStartScope(TRACE, "merge-points"); // iterate over all local blocks to give them unique ids. master.foreach ([&tolerance](ElementBlockT* b, // local block const diy::Master::ProxyWithLink&) -> void // communication proxy { b->MergeElements(tolerance); }); vtkLogEndScope("merge-points"); self->UpdateProgress(0.75); // now communicate point ownership information and assign ids to locally owned // points. vtkLogStartScope(TRACE, "exchange-ownership-ids"); diy::all_to_all(master, assigner, [](ElementBlockT* b, const diy::ReduceProxy& rp) { if (rp.round() == 0) { // now enqueue ownership information. b->EnqueueOwnershipInformation(rp); } else { // now dequeue ownership information and process locally to assign ids // to locally owned points and flag ghost points. b->DequeueOwnershipInformation(rp); } }); vtkLogEndScope("exchange-ownership-ids"); // exchange unique ids count so that we can determine global id offsets vtkLogStartScope(TRACE, "exchange-unique-ids"); diy::all_to_all(master, assigner, [](ElementBlockT* b, const diy::ReduceProxy& rp) { if (rp.round() == 0) { for (int i = rp.gid() + 1; i < rp.nblocks(); ++i) { rp.enqueue(rp.out_link().target(i), b->UniqueElementsCount); } } else { vtkIdType offset = 0; for (int src_gid = 0; src_gid < rp.gid(); ++src_gid) { vtkIdType msg; rp.dequeue(src_gid, msg); offset += msg; } b->AddOffset(offset); } }); vtkLogEndScope("exchange-unique-ids"); // exchange assigned ids. vtkLogStartScope(TRACE, "exchange-assigned-ids"); diy::all_to_all(master, assigner, [](ElementBlockT* b, const diy::ReduceProxy& rp) { if (rp.round() == 0) { b->EnqueueReplies(rp); } else { b->DequeueReplies(rp); } }); vtkLogEndScope("exchange-assigned-ids"); // final back communication to assign ids to ghosted points. vtkLogStartScope(TRACE, "exchange-ghosted-ids"); diy::all_to_all(master, assigner, [](ElementBlockT* b, const diy::ReduceProxy& rp) { if (rp.round() == 0) { b->EnqueueGhostedIds(rp); } else { b->DequeueGhostedIds(rp); } }); vtkLogEndScope("exchange-ghosted-ids"); self->UpdateProgress(1.0); return true; } VTK_ABI_NAMESPACE_END } namespace { struct CopyHiddenGhostPointsWorker { CopyHiddenGhostPointsWorker( vtkUnsignedCharArray* inputGhosts, vtkUnsignedCharArray* outputGhosts, bool cell_centers) : InputGhosts(inputGhosts) , OutputGhosts(outputGhosts) { if (cell_centers) { this->HiddenGhost = vtkDataSetAttributes::HIDDENCELL; } else { this->HiddenGhost = vtkDataSetAttributes::HIDDENPOINT; } } void operator()(vtkIdType startId, vtkIdType endId) { auto inputGhostsRange = vtk::DataArrayValueRange<1>(this->InputGhosts); auto outputGhostsRange = vtk::DataArrayValueRange<1>(this->OutputGhosts); for (vtkIdType id = startId; id < endId; ++id) { if (inputGhostsRange[id] & this->HiddenGhost) { outputGhostsRange[id] = inputGhostsRange[id]; } } } vtkUnsignedCharArray* InputGhosts; vtkUnsignedCharArray* OutputGhosts; unsigned char HiddenGhost; }; /** * This is the point type that keeps the coordinates for each point in the * dataset as well as enough information to track where that point came from so * that we can communicate back to the source once a unique global id has been * assigned. */ struct PointTT { static constexpr int attr_type = vtkDataObject::POINT; vtkTuple coords; int source_gid; vtkIdType source_id; // note: there's loss of precision here, but that's okay. this is only used by // DIY when building the kdtree float operator[](unsigned int index) const { return static_cast(this->coords[index]); } static std::vector GetElements(int gid, vtkPoints* pts, vtkDataSet*) { std::vector elems; elems.resize(pts->GetNumberOfPoints()); vtkSMPTools::For(0, pts->GetNumberOfPoints(), [&elems, pts, gid](vtkIdType start, vtkIdType end) { for (vtkIdType cc = start; cc < end; ++cc) { auto& pt = elems[cc]; pts->GetPoint(cc, pt.coords.GetData()); pt.source_gid = gid; pt.source_id = cc; } }); return elems; } static void Sort(std::vector& points) { // let's sort the points by source-id. This ensures that when a point is // duplicated among multiple blocks, the block with lower block-id owns the // point. Thus, keeping the numbering consistent. std::sort(points.begin(), points.end(), [](const PointTT& a, const PointTT& b) { return (a.source_gid == b.source_gid) ? (a.source_id < b.source_id) : (a.source_gid < b.source_gid); }); } static std::vector GenerateMergeMap( const std::vector& points, double tolerance) { std::vector mergemap(points.size(), -1); if (points.empty()) { return mergemap; } auto numPts = static_cast(points.size()); vtkNew grid; vtkNew pts; pts->SetDataTypeToDouble(); pts->SetNumberOfPoints(numPts); vtkSMPTools::For(0, numPts, [&](vtkIdType start, vtkIdType end) { for (vtkIdType cc = start; cc < end; ++cc) { pts->SetPoint(cc, points[cc].coords.GetData()); } }); grid->SetPoints(pts); vtkNew locator; locator->SetDataSet(grid); locator->SetTolerance(tolerance); locator->BuildLocator(); locator->MergePoints(tolerance, mergemap.data()); return mergemap; } }; struct CellTT { static constexpr int attr_type = vtkDataObject::CELL; vtkTuple center; int source_gid; vtkIdType source_id; std::vector point_ids; // note: there's loss of precision here, but that's okay. this is only used by // DIY when building the kdtree float operator[](unsigned int index) const { return static_cast(this->center[index]); } static std::vector GetElements(int gid, vtkPoints* centers, vtkDataSet* ds) { const vtkIdType ncells = ds->GetNumberOfCells(); // this is possible if vtkPoints created by vtkCellCenters included empty cells, which were // empty because they were ghosts. if (centers->GetNumberOfPoints() != ds->GetNumberOfCells()) { return std::vector(); } assert(centers->GetNumberOfPoints() == ncells); std::vector elems(ncells); vtkSMPThreadLocalObject tlIdList; if (ncells > 0) { // so that we can call GetCellPoints in vtkSMPTools::For ds->GetCellPoints(0, tlIdList.Local()); } auto pt_gids = vtkIdTypeArray::SafeDownCast(ds->GetPointData()->GetGlobalIds()); assert(ncells == 0 || pt_gids != nullptr); vtkSMPTools::For(0, ncells, [&](vtkIdType start, vtkIdType end) { auto ids = tlIdList.Local(); for (vtkIdType cc = start; cc < end; ++cc) { auto& cell = elems[cc]; centers->GetPoint(cc, cell.center.GetData()); cell.source_gid = gid; cell.source_id = cc; ds->GetCellPoints(cc, ids); cell.point_ids.resize(ids->GetNumberOfIds()); for (vtkIdType kk = 0, max = ids->GetNumberOfIds(); kk < max; ++kk) { cell.point_ids[kk] = pt_gids->GetTypedComponent(ids->GetId(kk), 0); } } }); return elems; } static void Sort(std::vector& cells) { // here, we are sorting such that for duplicated cells, we always order the // cell on the lower block before the one on the higher block. This is // essential to keep the cell numbering consistent. std::sort(cells.begin(), cells.end(), [](const CellTT& lhs, const CellTT& rhs) { return (lhs.point_ids == rhs.point_ids ? lhs.source_gid < rhs.source_gid : lhs.point_ids < rhs.point_ids); }); } static std::vector GenerateMergeMap( const std::vector& cells, double vtkNotUsed(tolerance)) { std::vector mergemap(cells.size(), -1); if (cells.empty()) { return mergemap; } mergemap[0] = 0; for (size_t cc = 1, max = cells.size(); cc < max; ++cc) { if (cells[cc - 1].point_ids == cells[cc].point_ids) { mergemap[cc] = mergemap[cc - 1]; } else { mergemap[cc] = static_cast(cc); } } return mergemap; } }; struct MessageItemTT { vtkIdType elem_id; vtkIdType index; }; template struct BlockT { private: void Enqueue(const diy::ReduceProxy& rp) { for (const auto& pair : this->OutMessage) { rp.enqueue(rp.out_link().target(pair.first), pair.second); } this->OutMessage.clear(); } public: vtkDataSet* Dataset{ nullptr }; std::vector Elements; std::vector MergeMap; vtkIdType UniqueElementsCount{ 0 }; std::map> OutMessage; vtkSmartPointer GlobalIds; vtkSmartPointer GhostArray; void Initialize(int self_gid, vtkPoints* points, vtkDataSet* dataset, bool cell_centers) { this->Dataset = dataset; this->Elements = ElementT::GetElements(self_gid, points, dataset); if (dataset) { unsigned char duplicateGhost = vtkDataSetAttributes::DUPLICATEPOINT; if (cell_centers) { duplicateGhost = vtkDataSetAttributes::DUPLICATECELL; } this->GlobalIds = vtkSmartPointer::New(); this->GlobalIds->SetName( ElementT::attr_type == vtkDataObject::POINT ? "GlobalPointIds" : "GlobalCellIds"); this->GlobalIds->SetNumberOfTuples(points->GetNumberOfPoints()); this->GlobalIds->FillValue(-1); dataset->GetAttributes(ElementT::attr_type)->SetGlobalIds(this->GlobalIds); this->GhostArray = vtkSmartPointer::New(); this->GhostArray->SetName(vtkDataSetAttributes::GhostArrayName()); this->GhostArray->SetNumberOfTuples(points->GetNumberOfPoints()); this->GhostArray->FillValue(duplicateGhost); if (vtkUnsignedCharArray* inputGhostPoints = dataset->GetAttributes(ElementT::attr_type)->GetGhostArray()) { ::CopyHiddenGhostPointsWorker worker(inputGhostPoints, this->GhostArray, cell_centers); vtkSMPTools::For(0, this->GhostArray->GetNumberOfValues(), worker); } // we're only adding ghost points, not cells. if (ElementT::attr_type == vtkDataObject::POINT) { dataset->GetAttributes(ElementT::attr_type)->AddArray(this->GhostArray); } } } void MergeElements(double tolerance) { // sort to make elements on lower gid's the primary elements ElementT::Sort(this->Elements); this->MergeMap = ElementT::GenerateMergeMap(this->Elements, tolerance); std::vector needs_replies(this->MergeMap.size()); for (size_t cc = 0, max = this->MergeMap.size(); cc < max; ++cc) { if (this->MergeMap[cc] != static_cast(cc)) { needs_replies[this->MergeMap[cc]] = 1; } } // populate out-message. for (size_t cc = 0, max = this->MergeMap.size(); cc < max; ++cc) { if (this->MergeMap[cc] == static_cast(cc)) { const auto& elem = this->Elements[cc]; MessageItemTT datum; datum.elem_id = elem.source_id; datum.index = needs_replies[cc] == 1 ? static_cast(cc) : vtkIdType(-1); this->OutMessage[elem.source_gid].push_back(datum); } } } void EnqueueOwnershipInformation(const diy::ReduceProxy& rp) { this->Enqueue(rp); } void DequeueOwnershipInformation(const diy::ReduceProxy& rp) { std::map> inmessage; for (int i = 0; i < rp.in_link().size(); ++i) { const int in_gid = rp.in_link().target(i).gid; while (rp.incoming(in_gid)) { auto& ownerships = inmessage[in_gid]; rp.dequeue(in_gid, ownerships); } } // we should not have received any message if we don't have a dataset. assert(this->Dataset != nullptr || inmessage.size() == 0); if (!this->Dataset) { return; } // We will tag any ghost that is not a hidden point in the input ghost array auto ghosts = vtk::DataArrayValueRange<1>(this->GhostArray); for (const auto& pair : inmessage) { for (const auto& data : pair.second) { if (!(ghosts[data.elem_id] & vtkDataSetAttributes::HIDDENPOINT)) { ghosts[data.elem_id] = 0; } } } // Assign global ids starting with 0 for locally owned elems. this->UniqueElementsCount = 0; for (vtkIdType cc = 0, max = this->GhostArray->GetNumberOfTuples(); cc < max; ++cc) { if (ghosts[cc] == 0) { this->GlobalIds->SetTypedComponent(cc, 0, this->UniqueElementsCount); this->UniqueElementsCount++; } } // Generate message send back assigned global ids to requesting blocks. for (const auto& pair : inmessage) { for (const auto& data : pair.second) { if (data.index != -1) { MessageItemTT reply; reply.index = data.index; reply.elem_id = this->GlobalIds->GetTypedComponent(data.elem_id, 0); this->OutMessage[pair.first].push_back(reply); } } } } void AddOffset(const vtkIdType offset) { if (this->GlobalIds == nullptr || offset == 0) { return; } vtkSMPTools::For(0, this->GlobalIds->GetNumberOfTuples(), [&offset, this](vtkIdType start, vtkIdType end) { for (vtkIdType cc = start; cc < end; ++cc) { const auto id = this->GlobalIds->GetTypedComponent(cc, 0); if (id != -1) { this->GlobalIds->SetTypedComponent(cc, 0, id + offset); } } }); // offset replies too. for (auto& pair : this->OutMessage) { for (auto& data : pair.second) { data.elem_id += offset; } } } void EnqueueReplies(const diy::ReduceProxy& rp) { this->Enqueue(rp); } void DequeueReplies(const diy::ReduceProxy& rp) { for (int i = 0; i < rp.in_link().size(); ++i) { const int in_gid = rp.in_link().target(i).gid; while (rp.incoming(in_gid)) { std::vector ownerships; rp.dequeue(in_gid, ownerships); for (const auto& data : ownerships) { // we're changing the id in our local storage to now be the global // id. this->Elements[data.index].source_id = data.elem_id; } } } const auto& mergemap = this->MergeMap; for (size_t cc = 0, max = mergemap.size(); cc < max; ++cc) { if (mergemap[cc] != static_cast(cc)) { auto& original_elem = this->Elements[mergemap[cc]]; auto& duplicate_elem = this->Elements[cc]; MessageItemTT data; data.elem_id = original_elem.source_id; data.index = duplicate_elem.source_id; this->OutMessage[duplicate_elem.source_gid].push_back(data); } } } void EnqueueGhostedIds(const diy::ReduceProxy& rp) { this->Enqueue(rp); } void DequeueGhostedIds(const diy::ReduceProxy& rp) { for (int i = 0; i < rp.in_link().size(); ++i) { const int in_gid = rp.in_link().target(i).gid; while (rp.incoming(in_gid)) { std::vector ownerships; rp.dequeue(in_gid, ownerships); assert(this->Dataset != nullptr || ownerships.size() == 0); for (const auto& data : ownerships) { this->GlobalIds->SetTypedComponent(data.index, 0, data.elem_id); } } } } }; } namespace diy { template <> struct Serialization<::CellTT> { static void save(BinaryBuffer& bb, const CellTT& c) { diy::save(bb, c.center); diy::save(bb, c.source_gid); diy::save(bb, c.source_id); diy::save(bb, c.point_ids); } static void load(BinaryBuffer& bb, CellTT& c) { c.point_ids.clear(); diy::load(bb, c.center); diy::load(bb, c.source_gid); diy::load(bb, c.source_id); diy::load(bb, c.point_ids); } }; } VTK_ABI_NAMESPACE_BEGIN vtkStandardNewMacro(vtkGenerateGlobalIds); vtkCxxSetObjectMacro(vtkGenerateGlobalIds, Controller, vtkMultiProcessController); //------------------------------------------------------------------------------ vtkGenerateGlobalIds::vtkGenerateGlobalIds() : Controller(nullptr) , Tolerance(0) { this->SetNumberOfInputPorts(1); this->SetNumberOfOutputPorts(1); this->SetController(vtkMultiProcessController::GetGlobalController()); } //------------------------------------------------------------------------------ vtkGenerateGlobalIds::~vtkGenerateGlobalIds() { this->SetController(nullptr); } //------------------------------------------------------------------------------ int vtkGenerateGlobalIds::RequestData( vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector) { auto inputDO = vtkDataObject::GetData(inputVector[0], 0); auto outputDO = vtkDataObject::GetData(outputVector, 0); outputDO->ShallowCopy(inputDO); // generate point ids first. { this->SetProgressShiftScale(0, 0.5); vtkLogScopeF(TRACE, "generate global point ids"); if (!impl::GenerateIds>(outputDO, this, false)) { this->SetProgressShiftScale(0, 1.0); return 0; } } // generate cell ids next (this needs global point ids) { this->SetProgressShiftScale(0.5, 0.5); vtkLogScopeF(TRACE, "generate global cell ids"); if (!impl::GenerateIds>(outputDO, this, true)) { this->SetProgressShiftScale(0, 1.0); return 0; } } this->SetProgressShiftScale(0, 1.0); return 1; } //------------------------------------------------------------------------------ void vtkGenerateGlobalIds::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os, indent); os << indent << "Controller: " << this->Controller << endl; os << indent << "Tolerance: " << this->Tolerance << endl; } VTK_ABI_NAMESPACE_END