Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Return fields required by TreeCompiler as vtk-m data set from filter.

Browse Source
This commit is contained in:
Gunther H. Weber 2020-09-21 19:05:10 -07:00
parent c4de6ba6e0
commit 7661b17e5b
6 changed files with 563 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
examples/contour_tree_distributed/CMakeLists.txt
View File

@ -59,7 +59,7 @@ find_package(VTKm REQUIRED QUIET)
# MPI
####################################
if (VTKm_ENABLE_MPI)
add_executable(ContourTree_Distributed ContourTreeApp.cxx)
add_executable(ContourTree_Distributed ContourTreeApp.cxx TreeCompiler.cxx)
target_link_libraries(ContourTree_Distributed vtkm_filter vtkm_io MPI::MPI_CXX)
vtkm_add_target_information(ContourTree_Distributed
MODIFY_CUDA_FLAGS
@ -70,6 +70,7 @@ if (VTKm_ENABLE_MPI)
mark_as_advanced(VTKM_EXAMPLE_CONTOURTREE_ENABLE_DEBUG_PRINT)
if (VTKM_EXAMPLE_CONTOURTREE_ENABLE_DEBUG_PRINT)
target_compile_definitions(ContourTree_Distributed PRIVATE "DEBUG_PRINT")
target_compile_definitions(ContourTree_Distributed PRIVATE "DEBUG_PRINT_CTUD")
endif()
if (TARGET vtkm::tbb)

89
examples/contour_tree_distributed/ContourTreeApp.cxx
View File

@ -98,8 +98,9 @@ VTKM_THIRDPARTY_POST_INCLUDE
#include <utility>
#include <vector>
#include "TreeCompiler.h"
#define PRESPLIT_FILE
using ValueType = vtkm::Float32;
namespace ctaug_ns = vtkm::worklet::contourtree_augmented;
@ -154,7 +155,42 @@ private:
std::vector<std::string> mCLOptions;
};
#if 0
struct PrintArrayContentsFunctor
{
template <typename T, typename S>
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<T, S>& array) const
{
this->PrintArrayPortal(array.GetPortalConstControl());
}
private:
template <typename PortalType >
VTKM_CONT void PrintArrayPortal(const PortalType& portal) const
{
for (vtkm::Id index = 0; index < portal.GetNumberOfValues(); index++)
{
// All ArrayPortal objects have ValueType for the type of each value.
using ValueType = typename PortalType::ValueType;
ValueType value = portal.Get(index);
vtkm::IdComponent numComponents =
vtkm::VecTraits <ValueType >::GetNumberOfComponents(value);
for (vtkm::IdComponent componentIndex = 0; componentIndex < numComponents; componentIndex ++)
{
std::cout << " " << vtkm::VecTraits<ValueType>::GetComponent(value, componentIndex);
}
//std::cout << std::endl;
}
}
};
template <typename VariantArrayType >
void PrintArrayContents(const VariantArrayType& array)
{
array.CastAndCall(PrintArrayContentsFunctor());
}
#endif
// Compute and render an isosurface for a uniform grid example
int main(int argc, char* argv[])
@ -471,6 +507,7 @@ int main(int argc, char* argv[])
}
// Read data
using ValueType = vtkm::FloatDefault;
std::vector<ValueType> values(numVertices);
if (filename.compare(filename.length() - 5, 5, ".bdem") == 0)
{
@ -534,6 +571,19 @@ int main(int argc, char* argv[])
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(nDims == 3 ? dims[2] : 0) });
#ifdef DEBUG_PRINT_CTUD
std::cout << "blockIndex: "
<< vtkm::Id3{ static_cast<vtkm::Id>(std::ceil(static_cast<float>(offset[0]) /
static_cast<float>(dims[0]))),
static_cast<vtkm::Id>(std::ceil(static_cast<float>(offset[1]) /
static_cast<float>(dims[1]))),
static_cast<vtkm::Id>(nDims == 3
? std::ceil(static_cast<float>(offset[2]) /
static_cast<float>(dims[2]))
: 0) }
<< std::endl;
#endif
if (blockNo == 0)
{
// FIXME: Hack: Approximate blocks per dim
@ -546,6 +596,9 @@ int main(int argc, char* argv[])
static_cast<float>(dims[2]))
: 1)
};
#ifdef DEBUG_PRINT_CTUD
std::cout << "blocksPerDim: " << blocksPerDim << std::endl;
#endif
}
}
@ -786,10 +839,40 @@ int main(int argc, char* argv[])
useMarchingCubes);
filter.SetActiveField("values");
// Execute the contour tree analysis. NOTE: If MPI is used the result will be
// a vtkm::cont::PartitionedDataSet instead of a vtkm::cont::DataSet
// Execute the contour tree analysis
auto result = filter.Execute(useDataSet);
#if 0
std::cout << "Result dataset has " << result.GetNumberOfPartitions() << " partitions" << std::endl;
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
auto ds = result.GetPartition(ds_no);
for (vtkm::Id f_no = 0; f_no < ds.GetNumberOfFields(); ++f_no)
{
auto field = ds.GetField(f_no);
std::cout << field.GetName() << ": ";
PrintArrayContents(field.GetData());
std::cout << std::endl;
}
}
#endif
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
TreeCompiler treeCompiler;
treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
char fname[256];
std::snprintf(fname,
sizeof(fname),
"TreeComplilerOutput_Rank%d_Block%d.dat",
rank,
static_cast<int>(ds_no));
FILE* out_file = std::fopen(fname, "wb");
treeCompiler.WriteBinary(out_file);
std::fclose(out_file);
}
currTime = totalTime.GetElapsedTime();
vtkm::Float64 computeContourTreeTime = currTime - prevTime;
prevTime = currTime;

239
examples/contour_tree_distributed/TreeCompiler.cxx Normal file
View File

@ -0,0 +1,239 @@
#include "TreeCompiler.h"
#include <iomanip>
#include <vtkm/cont/DataSet.h>
#define PRINT_WIDTH 12
// stream output
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node)
{ // stream output
outStream << node.lowEnd << " " << node.highEnd << " " << node.dataValue << " " << node.globalID
<< std::endl;
return outStream;
} // stream output
// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node)
{ // stream input
inStream >> node.lowEnd >> node.highEnd >> node.dataValue >> node.globalID;
return inStream;
} // stream input
// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
{ // TreeCompiler::AddHierarchicalTree()
// Copy relevant tree content to STL arrays
vtkm::cont::VariantArrayHandle dataValues_array = addedTree.GetField("DataValues").GetData();
std::vector<vtkm::FloatDefault> dataValues(dataValues_array.GetNumberOfValues());
auto dataValues_handle = vtkm::cont::make_ArrayHandle(dataValues, vtkm::CopyFlag::Off);
vtkm::cont::ArrayCopy(dataValues_array.ResetTypes(vtkm::List<vtkm::FloatDefault>{}),
dataValues_handle);
dataValues_handle.SyncControlArray();
auto regularNodeGlobalIds_array = addedTree.GetField("RegularNodeGlobalIds").GetData();
std::vector<vtkm::Id> regularNodeGlobalIds(regularNodeGlobalIds_array.GetNumberOfValues());
auto regularNodeGlobalIds_handle =
vtkm::cont::make_ArrayHandle(regularNodeGlobalIds, vtkm::CopyFlag::Off);
vtkm::cont::ArrayCopy(regularNodeGlobalIds_array.ResetTypes(vtkm::List<vtkm::Id>{}),
regularNodeGlobalIds_handle);
regularNodeGlobalIds_handle
.SyncControlArray(); //Forces values to get updated if copy happened on GPU
auto superarcs_array = addedTree.GetField("Superarcs").GetData();
std::vector<vtkm::Id> added_tree_superarcs(superarcs_array.GetNumberOfValues());
auto superarcs_handle = vtkm::cont::make_ArrayHandle(added_tree_superarcs, vtkm::CopyFlag::Off);
vtkm::cont::ArrayCopy(superarcs_array.ResetTypes(vtkm::List<vtkm::Id>{}), superarcs_handle);
superarcs_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU
auto supernodes_array = addedTree.GetField("Supernodes").GetData();
std::vector<vtkm::Id> added_tree_supernodes(supernodes_array.GetNumberOfValues());
auto supernodes_handle = vtkm::cont::make_ArrayHandle(added_tree_supernodes, vtkm::CopyFlag::Off);
vtkm::cont::ArrayCopy(supernodes_array.ResetTypes(vtkm::List<vtkm::Id>{}), supernodes_handle);
supernodes_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU
auto superparents_array = addedTree.GetField("Superparents").GetData();
std::vector<vtkm::Id> superparents(superparents_array.GetNumberOfValues());
auto superparents_handle = vtkm::cont::make_ArrayHandle(superparents, vtkm::CopyFlag::Off);
vtkm::cont::ArrayCopy(superparents_array.ResetTypes(vtkm::List<vtkm::Id>{}), superparents_handle);
superparents_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU
// loop through all of the supernodes in the hierarchical tree
for (indexType supernode = 0; supernode < added_tree_supernodes.size(); supernode++)
{ // per supernode
// retrieve the regular ID for the supernode
indexType regularId = added_tree_supernodes[supernode];
indexType globalId = regularNodeGlobalIds[regularId];
dataType dataVal = dataValues[regularId];
// retrieve the supernode at the far end
indexType superTo = added_tree_superarcs[supernode];
// now test - if it is NO_SUCH_ELEMENT, there are two possibilities
if (vtkm::worklet::contourtree_augmented::NoSuchElement(superTo))
{ // no Superto
// retrieve the superparent
indexType superparent = superparents[regularId];
// the root node will have itself as its superparent
if (superparent == supernode)
continue;
else
{ // not own superparent - an attachment point
// retrieve the superparent's from & to
indexType regularFrom = added_tree_supernodes[superparent];
indexType globalFrom = regularNodeGlobalIds[regularFrom];
indexType superParentTo = added_tree_superarcs[superparent];
indexType regularTo =
added_tree_supernodes[vtkm::worklet::contourtree_augmented::MaskedIndex(superParentTo)];
indexType globalTo = regularNodeGlobalIds[regularTo];
// test the superTo to see whether we ascend or descend
// note that we will never have NO_SUCH_ELEMENT here
if (vtkm::worklet::contourtree_augmented::IsAscending(superParentTo))
{ // ascending
this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalFrom, globalTo));
} // ascending
else
{ // descending
this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalFrom));
} // descending
} // not own superparent - an attachment point
} // no Superto
else
{ // there is a valid superarc
// retrieve the "to" and convert to global
indexType maskedTo = vtkm::worklet::contourtree_augmented::MaskedIndex(superTo);
indexType regularTo = added_tree_supernodes[maskedTo];
indexType globalTo = regularNodeGlobalIds[regularTo];
dataType dataTo = dataValues[regularTo];
// test the superTo to see whether we ascend or descend
// note that we will never have NO_SUCH_ELEMENT here
// we add both ends
if (vtkm::worklet::contourtree_augmented::IsAscending(superTo))
{ // ascending
this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalId, globalTo));
this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalId, globalTo));
} // ascending
else
{ // descending
this->supernodes.push_back(SupernodeOnSuperarc(globalId, dataVal, globalTo, globalId));
this->supernodes.push_back(SupernodeOnSuperarc(globalTo, dataTo, globalTo, globalId));
} // descending
} // there is a valid superarc
} // per supernode
} // TreeCompiler::AddHierarchicalTree()
// routine to compute the actual superarcs
void TreeCompiler::ComputeSuperarcs()
{ // TreeCompiler::ComputeSuperarcs()
// first we sort the vector
std::sort(supernodes.begin(), supernodes.end());
// we could do a unique test here, but it's easier just to suppress it inside the loop
// now we loop through it: note the -1
// this is because we know a priori that the last one is the last supernode on a superarc
// and would fail the test inside the loop. By putting it in the loop test, we avoid having
// to have an explicit if statement inside the loop
for (indexType supernode = 0; supernode < supernodes.size() - 1; supernode++)
{ // loop through supernodes
// this is actually painfully simple: if the (lowEnd, highEnd) don't match the next one,
// then we're at the end of the group and do nothing. Otherwise, we link to the next one
if ((supernodes[supernode].lowEnd != supernodes[supernode + 1].lowEnd) ||
(supernodes[supernode].highEnd != supernodes[supernode + 1].highEnd))
continue;
// if the supernode matches, then we have a repeat, and can suppress
if (supernodes[supernode].globalID == supernodes[supernode + 1].globalID)
continue;
// otherwise, add a superarc to the list
superarcs.push_back(Edge(supernodes[supernode].globalID, supernodes[supernode + 1].globalID));
} // loop through supernodes
// now sort them
std::sort(superarcs.begin(), superarcs.end());
} // TreeCompiler::ComputeSuperarcs()
// routine to print the superarcs
void TreeCompiler::PrintSuperarcs()
{ // TreeCompiler::PrintSuperarcs()
std::cout << "============" << std::endl;
std::cout << "Contour Tree" << std::endl;
for (indexType superarc = 0; superarc < superarcs.size(); superarc++)
{ // per superarc
if (superarcs[superarc].low < superarcs[superarc].high)
{ // order by ID not value
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << " ";
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << std::endl;
} // order by ID not value
else
{ // order by ID not value
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << " ";
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << std::endl;
} // order by ID not value
} // per superarc
} // TreeCompiler::PrintSuperarcs()
// routine to write out binary file
void TreeCompiler::WriteBinary(FILE* outFile)
{ // WriteBinary()
// do a bulk write of the entire contents
// no error checking, no type checking, no nothing
fwrite(&(supernodes[0]), sizeof(SupernodeOnSuperarc), supernodes.size(), outFile);
} // WriteBinary()
// routine to read in binary file and append
void TreeCompiler::ReadBinary(FILE* inFile)
{ // ReadBinary()
// use fseek to jump to the end
fseek(inFile, 0, SEEK_END);
// use fTell to retrieve the size of the file
long long nBytes = ftell(inFile);
// now rewind
rewind(inFile);
// compute how many elements are to be read
long long nSupernodes = nBytes / sizeof(SupernodeOnSuperarc);
// retrieve the current size
long long currentSize = supernodes.size();
// resize to add the right number
supernodes.resize(currentSize + nSupernodes);
// now read directly into the right chunk
fread(&(supernodes[currentSize]), sizeof(SupernodeOnSuperarc), nSupernodes, inFile);
} // ReadBinary()
// stream output - just dumps the supernodeonsuperarcs
std::ostream& operator<<(std::ostream& outStream, TreeCompiler& tree)
{ // stream output
for (indexType supernode = 0; supernode < tree.supernodes.size(); supernode++)
outStream << tree.supernodes[supernode];
return outStream;
} // stream output
// stream input - reads in the supernodeonsuperarcs & appends them
std::istream& operator>>(std::istream& inStream, TreeCompiler& tree)
{ // stream input
while (!inStream.eof())
{
SupernodeOnSuperarc tempNode;
inStream >> tempNode;
tree.supernodes.push_back(tempNode);
}
// we will overshoot, so subtract one
tree.supernodes.resize(tree.supernodes.size() - 1);
return inStream;
} // stream input

138
examples/contour_tree_distributed/TreeCompiler.h Normal file
View File

@ -0,0 +1,138 @@
#ifndef _TREECOMPILER_H_
#define _TREECOMPILER_H_
#include <iostream>
#include <vtkm/Types.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
// FIXME/HACK: Define here for compatibility with PPP TreeCompiler
typedef double dataType;
typedef unsigned long long indexType;
// small class for storing the contour arcs
class Edge
{ // Edge
public:
indexType low, high;
// constructor - defaults to -1
Edge(vtkm::Id Low = -1, vtkm::Id High = -1)
: low(Low)
, high(High)
{
}
}; // Edge
// comparison operator <
inline bool operator<(const Edge LHS, const Edge RHS)
{ // operator <
if (LHS.low < RHS.low)
return true;
if (LHS.low > RHS.low)
return false;
if (LHS.high < RHS.high)
return true;
if (LHS.high > RHS.high)
return false;
return false;
} // operator <
// a helper class which stores a single supernode inserted onto a superarc
class SupernodeOnSuperarc
{ // class SupernodeOnSuperarc
public:
// the global ID of the supernode
indexType globalID;
// the data value stored at the supernode
dataType dataValue;
// the low and high ends of the superarc it is on (may be itself)
indexType lowEnd, highEnd;
// constructor
SupernodeOnSuperarc(indexType GlobalID = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
dataType DataValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType LowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType HighEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)
: globalID(GlobalID)
, dataValue(DataValue)
, lowEnd(LowEnd)
, highEnd(HighEnd)
{ // constructor
} // constructor
}; // class SupernodeOnSuperarc
// overloaded comparison operator
// primary sort is by superarc (low, high),
// then secondary sort on datavalue
// tertiary on globalID to implement simulated simplicity
inline bool operator<(const SupernodeOnSuperarc& left, const SupernodeOnSuperarc& right)
{ // < operator
// simple lexicographic sort
if (left.lowEnd < right.lowEnd)
return true;
if (left.lowEnd > right.lowEnd)
return false;
if (left.highEnd < right.highEnd)
return true;
if (left.highEnd > right.highEnd)
return false;
if (left.dataValue < right.dataValue)
return true;
if (left.dataValue > right.dataValue)
return false;
if (left.globalID < right.globalID)
return true;
if (left.globalID > right.globalID)
return false;
// fall-through (shouldn't happen, but)
// if they're the same, it's false
return false;
} // < operator
// stream output
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node);
// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node);
// the class that compiles the contour tree
class TreeCompiler
{ // class TreeCompiler
public:
// we want a vector of supernodes on superarcs
std::vector<SupernodeOnSuperarc> supernodes;
// and a vector of Edges (the output)
std::vector<Edge> superarcs;
// default constructor sets it to empty
TreeCompiler()
{ // constructor
// clear out the supernode array
supernodes.resize(0);
// and the superarc array
superarcs.resize(0);
} // constructor
// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
void AddHierarchicalTree(const vtkm::cont::DataSet& addedTree);
// routine to compute the actual superarcs
void ComputeSuperarcs();
// routine to print the superarcs
void PrintSuperarcs();
// routine to write out binary file
void WriteBinary(FILE* outFile);
// routine to read in binary file & append to contents
void ReadBinary(FILE* inFile);
}; // class TreeCompiler
#endif

80
vtkm/filter/ContourTreeUniformDistributed.hxx
View File

@ -283,6 +283,7 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl(
boundaryTreeMaker.Construct(&localToGlobalIdRelabeler);
#ifdef DEBUG_PRINT_CTUD
// Print BRACT for debug
vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank();
char buffer[256];
@ -300,7 +301,7 @@ void ContourTreeUniformDistributed::ComputeLocalTreeImpl(
this->MultiBlockSpatialDecomposition.LocalBlockSizes.ReadPortal().Get(blockIndex),
this->MultiBlockSpatialDecomposition.GlobalSize)
<< std::endl;
#endif
// TODO: Fix calling conventions so dot print works here
// At this point, I'm reasonably certain that the contour tree has been computed regardless of data push/pull
// So although it might be logical to print things out earlier, I'll do it here
@ -378,7 +379,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::PreExecute(
template <typename FieldType, typename StorageType, typename DerivedPolicy>
VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
const vtkm::cont::PartitionedDataSet& input,
vtkm::cont::PartitionedDataSet&, // temporarily commented out: output,
vtkm::cont::PartitionedDataSet& result,
const vtkm::filter::FieldMetadata&, // dummy parameter for field meta data
const vtkm::cont::ArrayHandle<FieldType, StorageType>&, // dummy parameter to get the type
vtkm::filter::PolicyBase<DerivedPolicy>) // dummy parameter for policy
@ -440,6 +441,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
fieldData,
localGlobalMeshIndex);
#ifdef DEBUG_PRINT_CTUD
// ... save for debugging in text and .gv/.dot format
char buffer[256];
std::snprintf(buffer,
@ -460,6 +462,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
std::string(" Initial Step 5 BRACT Mesh"),
localDataBlocks[bi]->ContourTreeMeshes.back(),
worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL);
#endif
}
// Setup vtkmdiy to do global binary reduction of neighbouring blocks. See also RecuctionOperation struct for example
@ -533,9 +536,12 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
comm.barrier(); // Be safe!
std::vector<vtkm::cont::DataSet> hierarchicalTreeOutputDataSet(
localDataBlocks.size()); // DataSets for creating output data
master.foreach ([&](vtkm::worklet::contourtree_distributed::DistributedContourTreeBlockData<
FieldType>* b,
const vtkmdiy::Master::ProxyWithLink&) {
#ifdef DEBUG_PRINT_CTUD
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Fan In Complete" << std::endl
<< "# of CTs: " << b->ContourTrees.size() << std::endl
@ -564,6 +570,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
for (const auto& info : b->InteriorForests)
info.PrintContent(os);
os << std::endl;
#endif
// Fan out
auto nRounds = b->ContourTrees.size() - 1;
@ -619,6 +626,24 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
this->MultiBlockSpatialDecomposition.GlobalSize);
grafter.GraftInteriorForests(0, hierarchicalTree, fieldData, &localToGlobalIdRelabeler);
// Create data set from output
vtkm::cont::Field dataValuesField(
"DataValues", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.DataValues);
hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(dataValuesField);
vtkm::cont::Field regularNodeGlobalIdsField("RegularNodeGlobalIds",
vtkm::cont::Field::Association::WHOLE_MESH,
hierarchicalTree.RegularNodeGlobalIds);
hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(regularNodeGlobalIdsField);
vtkm::cont::Field superarcsField(
"Superarcs", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Superarcs);
hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(superarcsField);
vtkm::cont::Field supernodesField(
"Supernodes", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Supernodes);
hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(supernodesField);
vtkm::cont::Field superparentsField(
"Superparents", vtkm::cont::Field::Association::WHOLE_MESH, hierarchicalTree.Superparents);
hierarchicalTreeOutputDataSet[b->BlockIndex].AddField(superparentsField);
// TODO: GET THIS COMPILING
// save the corresponding .gv file
// std::string hierarchicalTreeFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(b->BlockIndex)) + "_Round_" + std::to_string(nRounds) + "_Hierarchical_Tree.gv";
@ -629,56 +654,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
// );
});
// FIXME: The following is still old code and needs to be changed!!!
#if 0
// Now run the contour tree algorithm on the last block to compute the final tree
vtkm::Id currNumIterations;
vtkm::worklet::contourtree_augmented::ContourTree currContourTree;
vtkm::worklet::contourtree_augmented::IdArrayType currSortOrder;
vtkm::worklet::ContourTreeAugmented worklet;
vtkm::cont::ArrayHandle<T> currField;
// Construct the mesh boundary exectuion object needed for boundary augmentation
vtkm::Id3 minIdx(0, 0, 0);
vtkm::Id3 maxIdx = this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize;
maxIdx[0] = maxIdx[0] - 1;
maxIdx[1] = maxIdx[1] - 1;
maxIdx[2] = maxIdx[2] > 0 ? (maxIdx[2] - 1) : 0;
// TODO Check setting. In parallel we should need to augment the tree in order to be able to do the merging
const unsigned int computeRegularStructure = 1;
auto meshBoundaryExecObj = localDataBlocks[0]->ContourTreeMeshes.back().GetMeshBoundaryExecutionObject(
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[0],
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[1],
minIdx,
maxIdx);
// Run the worklet to compute the final contour tree
worklet.Run(
localDataBlocks[0]->ContourTreeMeshes.back().SortedValues, // Unused param. Provide something to keep API happy
localDataBlocks[0]->ContourTreeMeshes.back(),
this->ContourTreeData,
this->MeshSortOrder,
currNumIterations,
computeRegularStructure,
meshBoundaryExecObj);
// Set the final mesh sort order we need to use
this->MeshSortOrder = localDataBlocks[0]->ContourTreeMeshes.back().GlobalMeshIndex;
// Remeber the number of iterations for the output
this->NumIterations = currNumIterations;
// Return the sorted values of the contour tree as the result
// TODO the result we return for the parallel and serial case are different right now. This should be made consistent. However, only in the parallel case are we useing the result output
vtkm::cont::DataSet temp;
vtkm::cont::Field rfield(this->GetOutputFieldName(),
vtkm::cont::Field::Association::WHOLE_MESH,
localDataBlocks[0]->ContourTreeMeshes.back().SortedValues);
temp.AddField(rfield);
output = vtkm::cont::PartitionedDataSet(temp);
/* Not necessary?
this->ContourTreeData = MultiBlockTreeHelper->LocalContourTrees[0];
this->MeshSortOrder = MultiBlockTreeHelper->LocalSortOrders[0];
*/
#endif
result = vtkm::cont::PartitionedDataSet(hierarchicalTreeOutputDataSet);
}

133
vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h
View File

@ -100,10 +100,11 @@ public:
) const
{
const auto selfid = rp.gid();
//DEBUGOUT//std::cout << "Entering ComputeDistributedContourTreeFunctor::operator() for GID " << selfid << std::endl;
#ifdef DEBUG_PRINT_CTUD
// Get rank (for debug output only)
const vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank();
#endif
// Here we do the deque first before the send due to the way the iteration is handled in DIY, i.e., in each iteration
// A block needs to first collect the data from its neighours and then send the combined block to its neighbours
@ -124,6 +125,9 @@ public:
rp.dequeue(ingid, otherBlockSize);
vtkm::worklet::contourtree_augmented::ContourTreeMesh<FieldType> otherContourTreeMesh;
rp.dequeue(ingid, otherContourTreeMesh);
#ifdef DEBUG_PRINT_CTUD
char buffer[256];
#if 0
// FIXME: Delete after debugging
std::cout << "Our block has extents: " << block->BlockOrigin << " " << block->BlockSize << std::endl;
@ -133,20 +137,19 @@ public:
block->ContourTreeMeshes.back().DebugPrint("OUR CTM", __FILE__, __LINE__);
otherContourTreeMesh.DebugPrint("OTHER CTM", __FILE__, __LINE__);
#endif
#if 0
char buffer[256];
std::snprintf(buffer, sizeof(buffer), "BeforeCombine_MyMesh_ContourTreeMesh_Rank%lld_Round%d.txt", rank, rp.round());
block->ContourTreeMeshes.back().Save(buffer);
std::snprintf(buffer, sizeof(buffer), "BeforeCombine_Other_ContourTreeMesh_Rank%lld_Round%d.txt", rank, rp.round());
otherContourTreeMesh.Save(buffer);
#endif
#endif
// Merge the two contour tree meshes
vtkm::cont::TryExecute(
MergeContourTreeMeshFunctor{}, otherContourTreeMesh, block->ContourTreeMeshes.back());
char buffer[256];
#ifdef DEBUG_PRINT_CTUD
std::snprintf(buffer,
sizeof(buffer),
"AfterCombine_ContourTreeMesh_Rank%d_Block%d_Round%d.txt",
@ -154,12 +157,12 @@ public:
static_cast<int>(block->BlockIndex),
rp.round());
block->ContourTreeMeshes.back().Save(buffer);
#if 0
// FIXME: Delete after debugging
std::cout << "================== AFTER COMBINING =================" << std::endl;
std::cout << "OUR CTM" << std::endl;
block->ContourTreeMeshes.back().DebugPrint("OUR CTM", __FILE__, __LINE__);
#endif
#endif
// Compute the origin and size of the new block
@ -196,13 +199,15 @@ public:
block->ContourTrees.back(),
currSortOrder,
currNumIterations,
1, // Fully augmented FIXME: Verify
1, // Fully augmented
meshBoundaryExecObj);
#ifdef DEBUG_PRINT_CTUD
#if 0
// FIXME: Delete after debugging
std::cout << "=================== BEGIN: COMBINED CONTOUR TREE =========================" << std::endl;
block->ContourTrees.back().PrintContent();
std::cout << "=================== END: COMBINED CONTOUR TREE =========================" << std::endl;
#endif
#endif
// Update block extents
@ -216,22 +221,19 @@ public:
// last round) then compute contour tree mesh to send and save it.
if (rp.round() != 0 && rp.out_link().size() != 0)
{
#ifdef DEBUG_PRINT_CTUD
char buffer[256];
#if 0
// If we have the fully augmented contour tree
block->ContourTreeMeshes.emplace_back(
block->ContourTrees.back().Arcs, block->ContourTreeMeshes.back()
);
char buffer[256];
std::snprintf(buffer, sizeof(buffer), "CombinedMeshes_GID%d_Round%d.txt", selfid, rp.round());
block->ContourTreeMeshes.back().Save(buffer);
#endif
#endif
vtkm::Id3 maxIdx{ block->BlockOrigin[0] + block->BlockSize[0] - 1,
block->BlockOrigin[1] + block->BlockSize[1] - 1,
block->BlockOrigin[2] + block->BlockSize[2] - 1 };
char buffer[256];
#ifdef DEBUG_PRINT_CTUD
std::snprintf(buffer,
sizeof(buffer),
"BRACTInputCTM_Rank%d_Block%d_Round%d.ctm_txt",
@ -256,28 +258,29 @@ public:
// TODO: GET THIS COMPILING
// save the corresponding .gv file for the contour tree mesh
// std::string contourTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_0_Combined_Mesh.gv");
// std::ofstream contourTreeMeshFile(contourTreeMeshFileName);
// contourTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint<FieldType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 0 Combined Mesh"),
// block->ContourTreeMeshes.back(), worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL);
//
// // and the ones for the contour tree regular and superstructures
// std::string regularStructureFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_1_Contour_Tree_Regular_Structure.gv");
// std::ofstream regularStructureFile(regularStructureFileName);
// regularStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint<T, MeshType, vtkm::worklet::contourtree_augmented::IdArrayType()>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 1 Contour Tree Regular Structure"),
// block->Meshes.back(),
// block->ContourTrees.back(),
// worklet::contourtree_distributed::SHOW_REGULAR_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS);
//
// std::string superStructureFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_2_Contour_Tree_Super_Structure.gv");
// std::ofstream superStructureFile(superStructureFileName);
// superStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint<T, MeshType, vtkm::worklet::contourtree_augmented::IdArrayType()>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 2 Contour Tree Super Structure"),
// block->Meshes.back(),
// block->ContourTrees.back(),
// worklet::contourtree_distributed::SHOW_SUPER_STRUCTURE|worklet::contourtree_distributed::SHOW_HYPER_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS|worklet::contourtree_distributed::SHOW_ALL_SUPERIDS|worklet::contourtree_distributed::SHOW_ALL_HYPERIDS);
// std::string contourTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_0_Combined_Mesh.gv");
// std::ofstream contourTreeMeshFile(contourTreeMeshFileName);
// contourTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint<FieldType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 0 Combined Mesh"),
// block->ContourTreeMeshes.back(), worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL);
//
// // and the ones for the contour tree regular and superstructures
// std::string regularStructureFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_1_Contour_Tree_Regular_Structure.gv");
// std::ofstream regularStructureFile(regularStructureFileName);
// regularStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint<T, MeshType, vtkm::worklet::contourtree_augmented::IdArrayType()>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 1 Contour Tree Regular Structure"),
// block->Meshes.back(),
// block->ContourTrees.back(),
// worklet::contourtree_distributed::SHOW_REGULAR_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS);
//
// std::string superStructureFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string("_Step_2_Contour_Tree_Super_Structure.gv");
// std::ofstream superStructureFile(superStructureFileName);
// superStructureFile << worklet::contourtree_distributed::ContourTreeDotGraphPrint<T, MeshType, vtkm::worklet::contourtree_augmented::IdArrayType()>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 2 Contour Tree Super Structure"),
// block->Meshes.back(),
// block->ContourTrees.back(),
// worklet::contourtree_distributed::SHOW_SUPER_STRUCTURE|worklet::contourtree_distributed::SHOW_HYPER_STRUCTURE|worklet::contourtree_distributed::SHOW_ALL_IDS|worklet::contourtree_distributed::SHOW_ALL_SUPERIDS|worklet::contourtree_distributed::SHOW_ALL_HYPERIDS);
#endif
// Compute BRACT
vtkm::worklet::contourtree_distributed::BoundaryTree boundaryTree;
@ -286,7 +289,6 @@ public:
this->GlobalSize, block->BlockOrigin, maxIdx);
// Make the BRACT and InteriorForest (i.e., residue)
block->InteriorForests.emplace_back();
std::cout << "Calling BRACT Maker" << std::endl;
auto boundaryTreeMaker = vtkm::worklet::contourtree_distributed::BoundaryTreeMaker<
vtkm::worklet::contourtree_augmented::ContourTreeMesh<FieldType>,
vtkm::worklet::contourtree_augmented::MeshBoundaryContourTreeMeshExec>(
@ -298,6 +300,7 @@ public:
// Construct the BRACT and InteriorForest. Since we are working on a ContourTreeMesh we do
// not need to provide and IdRelabeler here in order to compute the InteriorForest
boundaryTreeMaker.Construct();
#ifdef DEBUG_PRINT_CTUD
os << "================= OUTPUT ===================" << std::endl;
os << "+++++++++++++++++ BRACT +++++++++++++++++++" << std::endl;
os << boundaryTree.DebugPrint("validate", __FILE__, __LINE__);
@ -319,9 +322,12 @@ public:
rp.round());
std::ofstream dotStream(buffer);
dotStream << debug_dot;
#endif
// Construct contour tree mesh from BRACT
block->ContourTreeMeshes.emplace_back(
boundaryTree.VertexIndex, boundaryTree.Superarcs, block->ContourTreeMeshes.back());
#ifdef DEBUG_PRINT_CTUD
block->ContourTreeMeshes.back().PrintContent(os);
//char buffer[256];
std::snprintf(buffer,
@ -331,34 +337,35 @@ public:
static_cast<int>(block->BlockIndex),
block->ContourTreeMeshes.size() - 1);
block->ContourTreeMeshes.back().Save(buffer);
}
// TODO: GET THIS COMPILING
// save the Boundary Tree as a dot file
// std::string boundaryTreeFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_3_Boundary_Tree.gv");
// std::ofstream boundaryTreeFile(boundaryTreeFileName);
// boundaryTreeFile << vtkm::worklet::contourtree_distributed::BoundaryTreeDotGraphPrint
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 3 Boundary Tree"),
// block->Meshes.back()],
// block->BoundaryTrees.back());
//
// // and save the Interior Forest as another dot file
// std::string interiorForestFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_4_Interior_Forest.gv");
// std::ofstream interiorForestFile(interiorForestFileName);
// interiorForestFileName << InteriorForestDotGraphPrintFile<MeshType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 4 Interior Forest"),
// block->InteriorForests.back(),
// block->ContourTrees.back(),
// block->BoundaryTrees.back(),
// block->Meshes.back());
//
// // save the corresponding .gv file
// std::string boundaryTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_5_Boundary_Tree_Mesh.gv");
// std::ofstream boundaryTreeMeshFile(boundaryTreeMeshFileName);
// boundaryTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint<FieldType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 5 Boundary Tree Mesh"),
// block->ContourTreeMeshes.back(),
// worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL);
// TODO: GET THIS COMPILING
// save the Boundary Tree as a dot file
// std::string boundaryTreeFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_3_Boundary_Tree.gv");
// std::ofstream boundaryTreeFile(boundaryTreeFileName);
// boundaryTreeFile << vtkm::worklet::contourtree_distributed::BoundaryTreeDotGraphPrint
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 3 Boundary Tree"),
// block->Meshes.back()],
// block->BoundaryTrees.back());
//
// // and save the Interior Forest as another dot file
// std::string interiorForestFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_4_Interior_Forest.gv");
// std::ofstream interiorForestFile(interiorForestFileName);
// interiorForestFileName << InteriorForestDotGraphPrintFile<MeshType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 4 Interior Forest"),
// block->InteriorForests.back(),
// block->ContourTrees.back(),
// block->BoundaryTrees.back(),
// block->Meshes.back());
//
// // save the corresponding .gv file
// std::string boundaryTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast<int>(rank)) + std::string("_Block_") + std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_5_Boundary_Tree_Mesh.gv");
// std::ofstream boundaryTreeMeshFile(boundaryTreeMeshFileName);
// boundaryTreeMeshFile << vtkm::worklet::contourtree_distributed::ContourTreeMeshDotGraphPrint<FieldType>
// ( std::string("Block ") + std::to_string(static_cast<int>(block->BlockIndex)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 5 Boundary Tree Mesh"),
// block->ContourTreeMeshes.back(),
// worklet::contourtree_distributed::SHOW_CONTOUR_TREE_MESH_ALL);
#endif
}
// Send our current block (which is either our original block or the one we just combined from the ones we received) to our next neighbour.
// Once a rank has send his block (either in its orignal or merged form) it is done with the reduce