701 lines
25 KiB
C++
701 lines
25 KiB
C++
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: vtkDIYGhostUtilities.txx
|
|
|
|
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.
|
|
|
|
=========================================================================*/
|
|
#ifndef vtkDIYGhostUtilities_txx
|
|
#define vtkDIYGhostUtilities_txx
|
|
|
|
#include "vtkDIYGhostUtilities.h"
|
|
|
|
#include "vtkArrayDispatch.h"
|
|
#include "vtkBoundingBox.h"
|
|
#include "vtkCellData.h"
|
|
#include "vtkCommunicator.h"
|
|
#include "vtkDIYExplicitAssigner.h"
|
|
#include "vtkDIYUtilities.h"
|
|
#include "vtkDataArrayRange.h"
|
|
#include "vtkDataSetAttributes.h"
|
|
#include "vtkDataSetSurfaceFilter.h"
|
|
#include "vtkFieldData.h"
|
|
#include "vtkIdList.h"
|
|
#include "vtkImageData.h"
|
|
#include "vtkLogger.h"
|
|
#include "vtkMultiProcessController.h"
|
|
#include "vtkNew.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPoints.h"
|
|
#include "vtkPolyData.h"
|
|
#include "vtkRectilinearGrid.h"
|
|
#include "vtkSMPTools.h"
|
|
#include "vtkStructuredData.h"
|
|
#include "vtkStructuredGrid.h"
|
|
#include "vtkUnsignedCharArray.h"
|
|
#include "vtkUnstructuredGrid.h"
|
|
|
|
#include <atomic>
|
|
#include <functional>
|
|
#include <limits>
|
|
#include <string>
|
|
#include <type_traits>
|
|
#include <vector>
|
|
|
|
// clang-format off
|
|
#include "vtk_diy2.h"
|
|
#include VTK_DIY2(diy/assigner.hpp)
|
|
#include VTK_DIY2(diy/link.hpp)
|
|
#include VTK_DIY2(diy/master.hpp)
|
|
#include VTK_DIY2(diy/mpi.hpp)
|
|
#include VTK_DIY2(diy/partners/all-reduce.hpp)
|
|
#include VTK_DIY2(diy/reduce-operations.hpp)
|
|
// clang-format on
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<vtkImageData>
|
|
{
|
|
typedef ImageDataBlock BlockType;
|
|
static constexpr bool IsUnstructuredData = false;
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<vtkRectilinearGrid>
|
|
{
|
|
typedef RectilinearGridBlock BlockType;
|
|
static constexpr bool IsUnstructuredData = false;
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<vtkStructuredGrid>
|
|
{
|
|
typedef StructuredGridBlock BlockType;
|
|
static constexpr bool IsUnstructuredData = false;
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<vtkUnstructuredGrid>
|
|
{
|
|
typedef UnstructuredGridBlock BlockType;
|
|
static constexpr bool IsUnstructuredData = true;
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<vtkPolyData>
|
|
{
|
|
typedef PolyDataBlock BlockType;
|
|
static constexpr bool IsUnstructuredData = true;
|
|
};
|
|
|
|
namespace vtkDIYGhostUtilities_detail
|
|
{
|
|
//============================================================================
|
|
template <class ValueT, bool IsIntegerT = std::numeric_limits<ValueT>::is_integer>
|
|
struct Limits;
|
|
|
|
//============================================================================
|
|
template <class ValueT>
|
|
struct Limits<ValueT, true>
|
|
{
|
|
static constexpr ValueT Epsilon = 0;
|
|
static constexpr ValueT Min = 0;
|
|
};
|
|
|
|
//============================================================================
|
|
template <class ValueT>
|
|
struct Limits<ValueT, false>
|
|
{
|
|
static constexpr ValueT Epsilon = std::numeric_limits<ValueT>::epsilon();
|
|
static constexpr ValueT Min = std::numeric_limits<ValueT>::min();
|
|
};
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class ValueT>
|
|
ValueT ComputePrecision(ValueT val)
|
|
{
|
|
constexpr ValueT Epsilon = Limits<ValueT>::Epsilon;
|
|
constexpr ValueT Min = Limits<ValueT>::Min;
|
|
return std::max<ValueT>(val * Epsilon, Min);
|
|
}
|
|
|
|
//============================================================================
|
|
template <bool CellVersionT>
|
|
struct ConvertIdToCoordinatesFromCroppedExtentToFullExtentWorker;
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct ConvertIdToCoordinatesFromCroppedExtentToFullExtentWorker<true>
|
|
{
|
|
vtkIdType operator()(vtkIdType id, const vtkDIYGhostUtilities::ExtentType& fullExtent,
|
|
const vtkDIYGhostUtilities::ExtentType& croppedExtent)
|
|
{
|
|
int ijk[3];
|
|
vtkStructuredData::ComputeCellStructuredCoordsForExtent(id, croppedExtent.data(), ijk);
|
|
return vtkStructuredData::ComputeCellIdForExtent(fullExtent.data(), ijk);
|
|
}
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct ConvertIdToCoordinatesFromCroppedExtentToFullExtentWorker<false>
|
|
{
|
|
vtkIdType operator()(vtkIdType id, const vtkDIYGhostUtilities::ExtentType& fullExtent,
|
|
const vtkDIYGhostUtilities::ExtentType& croppedExtent)
|
|
{
|
|
int ijk[3];
|
|
vtkStructuredData::ComputePointStructuredCoordsForExtent(id, croppedExtent.data(), ijk);
|
|
return vtkStructuredData::ComputePointIdForExtent(fullExtent.data(), ijk);
|
|
}
|
|
};
|
|
|
|
//============================================================================
|
|
template <bool CellVersionT>
|
|
struct GhostFinder
|
|
{
|
|
using ExtentType = vtkDIYGhostUtilities::ExtentType;
|
|
|
|
GhostFinder(vtkUnsignedCharArray* ghosts, const ExtentType& fullExtent,
|
|
const ExtentType& croppedExtent, unsigned char ghostTarget)
|
|
: Ghosts(ghosts)
|
|
, FullExtent(fullExtent)
|
|
, CroppedExtent(croppedExtent)
|
|
, GhostTarget(ghostTarget)
|
|
, FoundGhost(false)
|
|
{
|
|
}
|
|
|
|
void operator()(vtkIdType startId, vtkIdType endId)
|
|
{
|
|
auto ghosts = vtk::DataArrayValueRange<1>(this->Ghosts);
|
|
using RefType = typename decltype(ghosts)::ReferenceType;
|
|
|
|
for (vtkIdType id = startId; id < endId; ++id)
|
|
{
|
|
if (this->FoundGhost)
|
|
{
|
|
return;
|
|
}
|
|
|
|
ConvertIdToCoordinatesFromCroppedExtentToFullExtentWorker<CellVersionT> worker;
|
|
vtkIdType remappedId = worker(id, this->FullExtent, this->CroppedExtent);
|
|
|
|
RefType ghost = ghosts[remappedId];
|
|
|
|
if (ghost & this->GhostTarget)
|
|
{
|
|
this->FoundGhost = true;
|
|
}
|
|
}
|
|
}
|
|
|
|
vtkUnsignedCharArray* Ghosts;
|
|
const ExtentType& FullExtent;
|
|
const ExtentType& CroppedExtent;
|
|
unsigned char GhostTarget;
|
|
|
|
std::atomic<bool> FoundGhost;
|
|
};
|
|
|
|
static constexpr unsigned char GHOST_CELL_BIT = 1;
|
|
static constexpr unsigned char GHOST_POINT_BIT = 2;
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class StructuredDataSetT>
|
|
unsigned char CleanGhostsReduceAllForStructuredData(
|
|
diy::Master& master, diy::Assigner& assigner, diy::RegularAllReducePartners& partners)
|
|
{
|
|
using BlockType =
|
|
typename vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<StructuredDataSetT>::BlockType;
|
|
using ExtentType = vtkDIYGhostUtilities::ExtentType;
|
|
|
|
unsigned char cleanGhostMask = 0;
|
|
|
|
diy::reduce(master, assigner, partners,
|
|
[&](BlockType* block, const diy::ReduceProxy& rp, const diy::RegularAllReducePartners&) {
|
|
// step 1 dequeue and merge
|
|
for (int i = 0; i < rp.in_link().size(); ++i)
|
|
{
|
|
int gid = rp.in_link().target(i).gid;
|
|
|
|
unsigned char currentMask;
|
|
rp.dequeue(gid, ¤tMask, 1);
|
|
cleanGhostMask |= currentMask;
|
|
}
|
|
|
|
unsigned char foundGhost = 0;
|
|
const ExtentType& extent = block->Information.Extent;
|
|
StructuredDataSetT* ds = block->Information.Input;
|
|
const int* dsExtent = ds->GetExtent();
|
|
const ExtentType fullExtent{ dsExtent[0], dsExtent[1], dsExtent[2], dsExtent[3], dsExtent[4],
|
|
dsExtent[5] };
|
|
if (vtkUnsignedCharArray* ghostCells = ds->GetCellGhostArray())
|
|
{
|
|
const vtkIdType numberOfCells = vtkStructuredData::GetNumberOfCells(extent.data());
|
|
vtkDIYGhostUtilities_detail::GhostFinder<true> worker(ghostCells, fullExtent, extent,
|
|
vtkDataSetAttributes::HIDDENCELL | vtkDataSetAttributes::DUPLICATECELL);
|
|
vtkSMPTools::For(0, numberOfCells, worker);
|
|
|
|
foundGhost |= GHOST_CELL_BIT * worker.FoundGhost;
|
|
}
|
|
if (vtkUnsignedCharArray* ghostPoints = ds->GetPointGhostArray())
|
|
{
|
|
const vtkIdType numberOfPoints = vtkStructuredData::GetNumberOfPoints(extent.data());
|
|
vtkDIYGhostUtilities_detail::GhostFinder<false> worker(ghostPoints, fullExtent, extent,
|
|
vtkDataSetAttributes::HIDDENPOINT | vtkDataSetAttributes::DUPLICATEPOINT);
|
|
vtkSMPTools::For(0, numberOfPoints, worker);
|
|
|
|
foundGhost |= GHOST_POINT_BIT * worker.FoundGhost;
|
|
}
|
|
|
|
// step 2 enqueue
|
|
for (int i = 0; i < rp.out_link().size(); ++i)
|
|
{
|
|
rp.enqueue(rp.out_link().target(i), &foundGhost, 1);
|
|
}
|
|
|
|
cleanGhostMask |= foundGhost;
|
|
});
|
|
|
|
return cleanGhostMask;
|
|
}
|
|
|
|
//============================================================================
|
|
template <bool IsUnstructuredDataT>
|
|
struct CleanGhostsReduceAllWorkerImpl;
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct CleanGhostsReduceAllWorkerImpl<true>
|
|
{
|
|
template <class DataSetT>
|
|
unsigned char operator()(diy::Master&, diy::Assigner&, diy::RegularAllReducePartners&)
|
|
{
|
|
return 0;
|
|
}
|
|
};
|
|
|
|
//============================================================================
|
|
template <>
|
|
struct CleanGhostsReduceAllWorkerImpl<false>
|
|
{
|
|
template <class DataSetT>
|
|
unsigned char operator()(
|
|
diy::Master& master, diy::Assigner& assigner, diy::RegularAllReducePartners& partners)
|
|
{
|
|
return vtkDIYGhostUtilities_detail::CleanGhostsReduceAllForStructuredData<DataSetT>(
|
|
master, assigner, partners);
|
|
}
|
|
};
|
|
|
|
//============================================================================
|
|
template <class DataSetT>
|
|
struct CleanGhostsReduceAllWorker
|
|
{
|
|
using DataSetInfo = vtkDIYGhostUtilities::DataSetTypeToBlockTypeConverter<DataSetT>;
|
|
using BlockType = typename DataSetInfo::BlockType;
|
|
static constexpr bool IsUnstructuredData = DataSetInfo::IsUnstructuredData;
|
|
|
|
unsigned char operator()(
|
|
diy::Master& master, diy::Assigner& assigner, diy::RegularAllReducePartners& partners)
|
|
{
|
|
return this->Impl.template operator()<DataSetT>(master, assigner, partners);
|
|
}
|
|
|
|
CleanGhostsReduceAllWorkerImpl<IsUnstructuredData> Impl;
|
|
};
|
|
|
|
//============================================================================
|
|
template <class DataSetT>
|
|
vtkSmartPointer<DataSetT> RemoveGhostArraysIfNeeded(
|
|
DataSetT* input, unsigned char ghostCleaningMask)
|
|
{
|
|
vtkSmartPointer<DataSetT> ds = vtkSmartPointer<DataSetT>::New();
|
|
ds->ShallowCopy(input);
|
|
if (!(ghostCleaningMask & GHOST_CELL_BIT))
|
|
{
|
|
ds->GetCellData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
|
|
}
|
|
if (!(ghostCleaningMask & GHOST_POINT_BIT))
|
|
{
|
|
ds->GetPointData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
|
|
}
|
|
return ds;
|
|
}
|
|
} // namesapce vtkDIYGhostUtilities_detail
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::ExchangeBoundingBoxes(
|
|
diy::Master& master, const vtkDIYExplicitAssigner& assigner, std::vector<DataSetT*>& inputs)
|
|
{
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
diy::all_to_all(
|
|
master, assigner, [&master, &inputs](BlockType* block, const diy::ReduceProxy& srp) {
|
|
int myBlockId = srp.gid();
|
|
int localId = master.lid(myBlockId);
|
|
auto& input = inputs[localId];
|
|
if (srp.round() == 0)
|
|
{
|
|
const double* bounds = input->GetBounds();
|
|
vtkBoundingBox& bb = block->BoundingBox;
|
|
bb = vtkBoundingBox(bounds);
|
|
|
|
vtkDIYGhostUtilities::InflateBoundingBoxIfNecessary(input, bb);
|
|
|
|
for (int i = 0; i < srp.out_link().size(); ++i)
|
|
{
|
|
const diy::BlockID& blockId = srp.out_link().target(i);
|
|
if (blockId.gid != myBlockId)
|
|
{
|
|
srp.enqueue(blockId, bb.GetMinPoint(), 3);
|
|
srp.enqueue(blockId, bb.GetMaxPoint(), 3);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
double minPoint[3], maxPoint[3];
|
|
for (int i = 0; i < static_cast<int>(srp.in_link().size()); ++i)
|
|
{
|
|
const diy::BlockID& blockId = srp.in_link().target(i);
|
|
if (blockId.gid != myBlockId)
|
|
{
|
|
srp.dequeue(blockId, minPoint, 3);
|
|
srp.dequeue(blockId, maxPoint, 3);
|
|
|
|
block->NeighborBoundingBoxes.emplace(blockId.gid,
|
|
vtkBoundingBox(
|
|
minPoint[0], maxPoint[0], minPoint[1], maxPoint[1], minPoint[2], maxPoint[2]));
|
|
}
|
|
}
|
|
}
|
|
});
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::ExchangeGhosts(diy::Master& master, std::vector<DataSetT*>& inputs)
|
|
{
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
master.foreach ([&master, &inputs](BlockType* block, const diy::Master::ProxyWithLink& cp) {
|
|
int myBlockId = cp.gid();
|
|
int localId = master.lid(myBlockId);
|
|
auto& input = inputs[localId];
|
|
|
|
for (int id = 0; id < static_cast<int>(cp.link()->size()); ++id)
|
|
{
|
|
const diy::BlockID& blockId = cp.link()->target(id);
|
|
vtkDIYGhostUtilities::EnqueueGhosts(cp, blockId, input, block);
|
|
}
|
|
});
|
|
|
|
master.exchange();
|
|
|
|
master.foreach ([](BlockType* block, const diy::Master::ProxyWithLink& cp) {
|
|
std::vector<int> incoming;
|
|
cp.incoming(incoming);
|
|
for (const int& gid : incoming)
|
|
{
|
|
// we need this extra check because incoming is not empty when using only one block
|
|
if (!cp.incoming(gid).empty())
|
|
{
|
|
vtkDIYGhostUtilities::DequeueGhosts(cp, gid, block->BlockStructures.at(gid));
|
|
}
|
|
}
|
|
});
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class BlockT>
|
|
vtkDIYGhostUtilities::LinkMap vtkDIYGhostUtilities::ComputeLinkMapUsingBoundingBoxes(
|
|
const diy::Master& master)
|
|
{
|
|
LinkMap linkMap(master.size());
|
|
|
|
for (int localId = 0; localId < static_cast<int>(master.size()); ++localId)
|
|
{
|
|
Links& links = linkMap[localId];
|
|
BlockT* block = master.block<BlockT>(localId);
|
|
vtkBoundingBox& localbb = block->BoundingBox;
|
|
BlockMapType<vtkBoundingBox>& bb = block->NeighborBoundingBoxes;
|
|
for (auto item : bb)
|
|
{
|
|
if (localbb.Intersects(item.second))
|
|
{
|
|
links.insert(item.first);
|
|
}
|
|
}
|
|
}
|
|
return linkMap;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::CopyInputsAndAllocateGhosts(diy::Master& master, diy::Assigner& assigner,
|
|
diy::RegularAllReducePartners& partners, std::vector<DataSetT*>& inputs,
|
|
std::vector<DataSetT*>& outputs, int outputGhostLevels)
|
|
{
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
vtkDIYGhostUtilities_detail::CleanGhostsReduceAllWorker<DataSetT> cleaner;
|
|
unsigned char ghostCleaningMask = cleaner(master, assigner, partners);
|
|
|
|
for (int localId = 0; localId < static_cast<int>(inputs.size()); ++localId)
|
|
{
|
|
DataSetT* input = inputs[localId];
|
|
DataSetT* output = outputs[localId];
|
|
|
|
BlockType* block = master.block<BlockType>(localId);
|
|
vtkSmartPointer<DataSetT> cleanedInput =
|
|
vtkDIYGhostUtilities_detail::RemoveGhostArraysIfNeeded(input, ghostCleaningMask);
|
|
|
|
// If we are isolated (no connection with other blocks), just shallow-copy
|
|
if (block->BlockStructures.empty())
|
|
{
|
|
output->ShallowCopy(cleanedInput);
|
|
}
|
|
// If we fetch 0 levels of ghosts AND there were no ghost cells in the input,
|
|
// we can partially shallow copy (we can't shallow copy points as the interfaces
|
|
// can be written over by other blocks)
|
|
else if (outputGhostLevels == 0 && !block->Information.InputNeedsGhostsPeeledOff() &&
|
|
!(ghostCleaningMask & vtkDIYGhostUtilities_detail::GHOST_CELL_BIT))
|
|
{
|
|
output->CopyStructure(input);
|
|
output->GetPointData()->DeepCopy(cleanedInput->GetPointData());
|
|
output->GetCellData()->ShallowCopy(cleanedInput->GetCellData());
|
|
output->GetFieldData()->ShallowCopy(input->GetFieldData());
|
|
}
|
|
// In the general case, deep copy the input and allocate the geometry
|
|
// for the new ghost cells.
|
|
else
|
|
{
|
|
vtkDIYGhostUtilities::DeepCopyInputAndAllocateGhosts(block, cleanedInput, output);
|
|
}
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::AddGhostArrays(diy::Master& master, std::vector<DataSetT*>& outputs)
|
|
{
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
for (int localId = 0; localId < static_cast<int>(outputs.size()); ++localId)
|
|
{
|
|
DataSetT* output = outputs[localId];
|
|
BlockType* block = master.block<BlockType>(localId);
|
|
|
|
output->GetPointData()->AddArray(block->GhostPointArray);
|
|
|
|
if (block->GhostCellArray)
|
|
{
|
|
output->GetCellData()->AddArray(block->GhostCellArray);
|
|
}
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::InitializeGhostArrays(
|
|
diy::Master& master, std::vector<DataSetT*>& outputs, int outputGhostLevels)
|
|
{
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
for (int localId = 0; localId < static_cast<int>(outputs.size()); ++localId)
|
|
{
|
|
DataSetT* output = outputs[localId];
|
|
BlockType* block = master.block<BlockType>(localId);
|
|
|
|
if (outputGhostLevels != 0)
|
|
{
|
|
vtkDIYGhostUtilities::InitializeGhostCellArray(block, output);
|
|
}
|
|
|
|
vtkDIYGhostUtilities::InitializeGhostPointArray(block, output);
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::InitializeGhostCellArray(
|
|
typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType* block, DataSetT* output)
|
|
{
|
|
vtkSmartPointer<vtkUnsignedCharArray>& ghostCellArray = block->GhostCellArray;
|
|
ghostCellArray = output->GetGhostArray(vtkDataObject::FIELD_ASSOCIATION_CELLS);
|
|
if (!ghostCellArray)
|
|
{
|
|
ghostCellArray = vtkSmartPointer<vtkUnsignedCharArray>::New();
|
|
ghostCellArray->SetName(vtkDataSetAttributes::GhostArrayName());
|
|
ghostCellArray->SetNumberOfComponents(1);
|
|
ghostCellArray->SetNumberOfValues(output->GetNumberOfCells());
|
|
ghostCellArray->Fill(0);
|
|
}
|
|
else
|
|
{
|
|
vtkDIYGhostUtilities::ReinitializeSelectedBits(
|
|
ghostCellArray, vtkDataSetAttributes::DUPLICATECELL);
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class DataSetT>
|
|
void vtkDIYGhostUtilities::InitializeGhostPointArray(
|
|
typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType* block, DataSetT* output)
|
|
{
|
|
vtkSmartPointer<vtkUnsignedCharArray>& ghostPointArray = block->GhostPointArray;
|
|
ghostPointArray = output->GetGhostArray(vtkDataObject::FIELD_ASSOCIATION_POINTS);
|
|
if (!ghostPointArray)
|
|
{
|
|
ghostPointArray = vtkSmartPointer<vtkUnsignedCharArray>::New();
|
|
ghostPointArray->SetName(vtkDataSetAttributes::GhostArrayName());
|
|
ghostPointArray->SetNumberOfComponents(1);
|
|
ghostPointArray->SetNumberOfValues(output->GetNumberOfPoints());
|
|
ghostPointArray->Fill(0);
|
|
}
|
|
else
|
|
{
|
|
vtkDIYGhostUtilities::ReinitializeSelectedBits(
|
|
ghostPointArray, vtkDataSetAttributes::DUPLICATEPOINT);
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
/**
|
|
* Main pipeline generating ghosts.
|
|
*/
|
|
template <class DataSetT>
|
|
int vtkDIYGhostUtilities::GenerateGhostCells(std::vector<DataSetT*>& inputs,
|
|
std::vector<DataSetT*>& outputs, int outputGhostLevels, vtkMultiProcessController* controller)
|
|
{
|
|
static_assert((std::is_base_of<vtkImageData, DataSetT>::value ||
|
|
std::is_base_of<vtkRectilinearGrid, DataSetT>::value ||
|
|
std::is_base_of<vtkStructuredGrid, DataSetT>::value ||
|
|
std::is_base_of<vtkUnstructuredGrid, DataSetT>::value ||
|
|
std::is_base_of<vtkPolyData, DataSetT>::value),
|
|
"Input data set type is not supported.");
|
|
|
|
using BlockType = typename DataSetTypeToBlockTypeConverter<DataSetT>::BlockType;
|
|
|
|
const int size = static_cast<int>(inputs.size());
|
|
if (size != static_cast<int>(outputs.size()))
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
std::string logMessage = size
|
|
? std::string("Generating ghosts for ") + std::string(outputs[0]->GetClassName())
|
|
: std::string("No ghosts to generate for empty rank");
|
|
vtkLogStartScope(TRACE, logMessage.c_str());
|
|
|
|
vtkDIYGhostUtilities::CloneGeometricStructures(inputs, outputs);
|
|
|
|
vtkLogStartScope(TRACE, "Instantiating diy communicator");
|
|
diy::mpi::communicator comm = vtkDIYUtilities::GetCommunicator(controller);
|
|
vtkLogEndScope("Instantiating diy communicator");
|
|
|
|
vtkLogStartScope(TRACE, "Instantiating master");
|
|
diy::Master master(
|
|
comm, 1, -1, []() { return static_cast<void*>(new BlockType()); },
|
|
[](void* b) -> void { delete static_cast<BlockType*>(b); });
|
|
vtkLogEndScope("Instantiating master");
|
|
|
|
vtkLogStartScope(TRACE, "Instantiating assigner");
|
|
vtkDIYExplicitAssigner assigner(comm, size);
|
|
vtkLogEndScope("Instantiating assigner");
|
|
|
|
if (!size)
|
|
{
|
|
// In such instance, we can just terminate. We are empty an finished communicating with other
|
|
// ranks.
|
|
vtkLogEndScope(logMessage.c_str());
|
|
return 1;
|
|
}
|
|
|
|
vtkLogStartScope(TRACE, "Decomposing master");
|
|
diy::RegularDecomposer<diy::DiscreteBounds> decomposer(
|
|
/*dim*/ 1, diy::interval(0, assigner.nblocks() - 1), assigner.nblocks());
|
|
decomposer.decompose(comm.rank(), assigner, master);
|
|
vtkLogEndScope("Decomposing master");
|
|
|
|
diy::RegularAllReducePartners partners(decomposer, 2);
|
|
|
|
// At this step, we gather data from the inputs and store it inside the local blocks
|
|
// so we don't have to carry extra parameters later.
|
|
vtkLogStartScope(TRACE, "Setup block self information.");
|
|
vtkDIYGhostUtilities::InitializeBlocks(master, inputs);
|
|
vtkLogEndScope("Setup block self information.");
|
|
|
|
vtkLogStartScope(TRACE, "Exchanging bounding boxes");
|
|
vtkDIYGhostUtilities::ExchangeBoundingBoxes(master, assigner, inputs);
|
|
vtkLogEndScope("Exchanging bounding boxes");
|
|
|
|
// We compute a temporary link map that weeds out data sets that do not have
|
|
// overlapping bounding boxes.
|
|
vtkLogStartScope(TRACE, "Computing temporary link map using bounding boxes.");
|
|
LinkMap temporaryLinkMap =
|
|
vtkDIYGhostUtilities::ComputeLinkMapUsingBoundingBoxes<BlockType>(master);
|
|
vtkLogEndScope("Computing temporary link map using bounding boxes.");
|
|
|
|
// We link our blocks using the temporary link map
|
|
vtkLogStartScope(TRACE, "Relinking blocks using temporary link map");
|
|
vtkDIYUtilities::Link(master, assigner, temporaryLinkMap);
|
|
vtkLogEndScope("Relinking blocks using temporary link map");
|
|
|
|
// Here, we exchange structural information between blocks that will be used to
|
|
// determine if blocks are actually adjacent or not.
|
|
vtkLogStartScope(TRACE, "Exchanging block structures");
|
|
vtkDIYGhostUtilities::ExchangeBlockStructures(master, inputs);
|
|
vtkLogEndScope("Exchanging block structures");
|
|
|
|
// The structural information that has been exchanged is used to compute
|
|
// the final link map, mapping blocks that will actually exchange ghosts.
|
|
vtkLogStartScope(TRACE, "Creating link map between connected blocks");
|
|
LinkMap linkMap = vtkDIYGhostUtilities::ComputeLinkMap(master, inputs, outputGhostLevels);
|
|
vtkLogEndScope("Creating link map between connected blocks");
|
|
|
|
vtkLogStartScope(TRACE, "Relinking blocks using link map");
|
|
vtkDIYUtilities::Link(master, assigner, linkMap);
|
|
vtkLogEndScope("Relinking blocks using link map");
|
|
|
|
vtkLogStartScope(TRACE, "Exchanging ghost data between blocks");
|
|
vtkDIYGhostUtilities::ExchangeGhosts(master, inputs);
|
|
vtkLogEndScope("Exchanging ghost data between blocks");
|
|
|
|
vtkLogStartScope(TRACE, "Allocating ghosts in outputs");
|
|
vtkDIYGhostUtilities::CopyInputsAndAllocateGhosts(
|
|
master, assigner, partners, inputs, outputs, outputGhostLevels);
|
|
vtkLogEndScope("Allocating ghosts in outputs");
|
|
|
|
vtkLogStartScope(TRACE, "Initializing ghost arrays in outputs");
|
|
vtkDIYGhostUtilities::InitializeGhostArrays(master, outputs, outputGhostLevels);
|
|
vtkLogEndScope("Initializing ghost arrays in outputs");
|
|
|
|
vtkLogStartScope(TRACE, "Filling local ghosts with received data from other blocks");
|
|
vtkDIYGhostUtilities::FillGhostArrays(master, outputs, outputGhostLevels);
|
|
vtkLogEndScope("Filling local ghosts with received data from other blocks");
|
|
|
|
vtkLogStartScope(TRACE, "Adding ghost arrays to point and / or cell data");
|
|
vtkDIYGhostUtilities::AddGhostArrays(master, outputs);
|
|
vtkLogEndScope("Adding ghost arrays to point and / or cell data");
|
|
|
|
vtkLogEndScope(logMessage.c_str());
|
|
|
|
return 1;
|
|
}
|
|
|
|
#endif
|