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

Refactor mesh classes: switch from row/col/slice to vtkm::Id2/Id3

Browse Source
This commit is contained in:
Gunther H. Weber 2020-09-15 09:11:00 -07:00
parent d00071b4e1
commit 991f7a85ad
55 changed files with 1329 additions and 1004 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
examples/contour_tree_augmented/ContourTreeApp.cxx
View File

@ -387,8 +387,8 @@ int main(int argc, char* argv[])
#ifdef WITH_MPI
#ifdef DEBUG_PRINT
// From https://www.unix.com/302983597-post2.html
char* cstr_filename = new char[15];
snprintf(cstr_filename, sizeof(filename), "cout_%d.log", rank);
char cstr_filename[32];
snprintf(cstr_filename, sizeof(cstr_filename), "cout_%d.log", rank);
int out = open(cstr_filename, O_RDWR | O_CREAT | O_APPEND, 0600);
if (-1 == out)
{
@ -417,8 +417,6 @@ int main(int argc, char* argv[])
perror("cannot redirect stderr");
return 255;
}
delete[] cstr_filename;
#endif
#endif
@ -444,12 +442,12 @@ int main(int argc, char* argv[])
// Copy the data into the values array so we can construct a multiblock dataset
// TODO All we should need to do to implement BOV support is to copy the values
// in the values vector and copy the dimensions in the dims vector
vtkm::Id nRows, nCols, nSlices;
vtkm::worklet::contourtree_augmented::GetRowsColsSlices temp;
temp(inDataSet.GetCellSet(), nRows, nCols, nSlices);
dims[0] = nRows;
dims[1] = nCols;
dims[2] = nSlices;
vtkm::Id3 meshSize;
vtkm::worklet::contourtree_augmented::GetPointDimensions temp;
temp(inDataSet.GetCellSet(), meshSize);
dims[0] = meshSize[0];
dims[1] = meshSize[1];
dims[2] = meshSize[2];
auto tempField = inDataSet.GetField("values").GetData();
values.resize(static_cast<std::size_t>(tempField.GetNumberOfValues()));
auto tempFieldHandle = tempField.AsVirtual<ValueType>().ReadPortal();
@ -515,6 +513,13 @@ int main(int argc, char* argv[])
dataReadTime = currTime - prevTime;
prevTime = currTime;
// swap dims order
{
vtkm::Id temp = dims[0];
dims[0] = dims[1];
dims[1] = temp;
}
#ifndef WITH_MPI // We only need the inDataSet if are not using MPI otherwise we'll constructe a multi-block dataset
// build the input dataset
vtkm::cont::DataSetBuilderUniform dsb;
@ -522,16 +527,16 @@ int main(int argc, char* argv[])
if (nDims == 2)
{
vtkm::Id2 vdims;
vdims[0] = static_cast<vtkm::Id>(dims[1]);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(dims[1]);
inDataSet = dsb.Create(vdims);
}
// 3D data
else
{
vtkm::Id3 vdims;
vdims[0] = static_cast<vtkm::Id>(dims[1]);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(dims[1]);
vdims[2] = static_cast<vtkm::Id>(dims[2]);
inDataSet = dsb.Create(vdims);
}
@ -594,8 +599,8 @@ int main(int argc, char* argv[])
{
VTKM_LOG_IF_S(vtkm::cont::LogLevel::Error,
rank == 0,
"Number of ranks to large for data. Use " << lastDimSize / 2
<< "or fewer ranks");
"Number of ranks too large for data. Use " << lastDimSize / 2
<< "or fewer ranks");
MPI_Finalize();
return EXIT_FAILURE;
}
@ -629,8 +634,8 @@ int main(int argc, char* argv[])
if (nDims == 2)
{
vtkm::Id2 vdims;
vdims[0] = static_cast<vtkm::Id>(currBlockSize);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(currBlockSize);
vtkm::Vec<ValueType, 2> origin(0, blockIndex * blockSize);
vtkm::Vec<ValueType, 2> spacing(1, 1);
ds = dsb.Create(vdims, origin, spacing);
@ -645,8 +650,8 @@ int main(int argc, char* argv[])
else
{
vtkm::Id3 vdims;
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(dims[1]);
vdims[0] = static_cast<vtkm::Id>(dims[1]);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[2] = static_cast<vtkm::Id>(currBlockSize);
vtkm::Vec<ValueType, 3> origin(0, 0, (blockIndex * blockSize));
vtkm::Vec<ValueType, 3> spacing(1, 1, 1);
@ -690,6 +695,21 @@ int main(int argc, char* argv[])
vtkm::Float64 computeContourTreeTime = currTime - prevTime;
prevTime = currTime;
#ifdef WITH_MPI
#ifdef DEBUG_PRINT
std::cout << std::flush;
close(out);
std::cerr << std::flush;
close(err);
dup2(save_out, fileno(stdout));
dup2(save_err, fileno(stderr));
close(save_out);
close(save_err);
#endif
#endif
////////////////////////////////////////////
// Compute the branch decomposition
////////////////////////////////////////////

62
vtkm/cont/testing/MakeTestDataSet.h
View File

@ -43,12 +43,14 @@ public:
vtkm::cont::DataSet Make2DUniformDataSet0();
vtkm::cont::DataSet Make2DUniformDataSet1();
vtkm::cont::DataSet Make2DUniformDataSet2();
vtkm::cont::DataSet Make2DUniformDataSet3();
// 3D uniform datasets.
vtkm::cont::DataSet Make3DUniformDataSet0();
vtkm::cont::DataSet Make3DUniformDataSet1();
vtkm::cont::DataSet Make3DUniformDataSet2();
vtkm::cont::DataSet Make3DUniformDataSet3(const vtkm::Id3 dims = vtkm::Id3(10));
vtkm::cont::DataSet Make3DUniformDataSet4();
vtkm::cont::DataSet Make3DRegularDataSet0();
vtkm::cont::DataSet Make3DRegularDataSet1();
@ -177,7 +179,7 @@ inline vtkm::cont::DataSet MakeTestDataSet::Make2DUniformDataSet0()
return dataSet;
}
//Make a simple 2D, 16 cell uniform dataset.
//Make a simple 2D, 16 cell uniform dataset (5x5.txt)
inline vtkm::cont::DataSet MakeTestDataSet::Make2DUniformDataSet1()
{
vtkm::cont::DataSetBuilderUniform dsb;
@ -233,6 +235,29 @@ inline vtkm::cont::DataSet MakeTestDataSet::Make2DUniformDataSet2()
return dataSet;
}
//Make a simple 2D, 56 cell uniform dataset. (8x9test.txt)
inline vtkm::cont::DataSet MakeTestDataSet::Make2DUniformDataSet3()
{
vtkm::cont::DataSetBuilderUniform dsb;
constexpr vtkm::Id2 dimensions(9, 8);
vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
constexpr vtkm::Id nVerts = 72;
constexpr vtkm::Float32 pointvar[nVerts] = {
29.0f, 37.0f, 39.0f, 70.0f, 74.0f, 84.0f, 38.0f, 36.0f, 26.0f, 27.0f, 100.0f, 49.0f,
72.0f, 85.0f, 89.0f, 83.0f, 28.0f, 24.0f, 25.0f, 47.0f, 50.0f, 73.0f, 86.0f, 90.0f,
71.0f, 82.0f, 22.0f, 23.0f, 75.0f, 79.0f, 48.0f, 69.0f, 87.0f, 88.0f, 81.0f, 18.0f,
19.0f, 76.0f, 80.0f, 78.0f, 46.0f, 68.0f, 67.0f, 40.0f, 16.0f, 17.0f, 41.0f, 77.0f,
45.0f, 35.0f, 20.0f, 21.0f, 32.0f, 15.0f, 13.0f, 42.0f, 43.0f, 44.0f, 34.0f, 33.0f,
31.0f, 30.0f, 14.0f, 12.0f, 11.0f, 10.0f, 9.0f, 8.0f, 7.0f, 6.0f, 5.0f, 0.0f
};
dataSet.AddPointField("pointvar", pointvar, nVerts);
return dataSet;
}
//Make a simple 3D, 4 cell uniform dataset.
inline vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet0()
{
@ -254,7 +279,7 @@ inline vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet0()
return dataSet;
}
//Make a simple 3D, 64 cell uniform dataset.
//Make a simple 3D, 64 cell uniform dataset. (5b 5x5x5)
inline vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet1()
{
vtkm::cont::DataSetBuilderUniform dsb;
@ -372,6 +397,39 @@ inline vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet3(const vtkm::Id
return dataSet;
}
//Make a simple 3D, 120 cell uniform dataset. (This is the data set from
//Make3DUniformDataSet1 upsampled from 5x5x to 5x6x7.)
inline vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet4()
{
vtkm::cont::DataSetBuilderUniform dsb;
constexpr vtkm::Id3 dimensions(5, 6, 7);
vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
constexpr vtkm::Id nVerts = 210;
constexpr vtkm::Float32 pointvar[nVerts] = {
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.53f, 0.48f, 0.45f,
0.0f, 0.0f, 0.64f, 0.56f, 0.61f, 0.0f, 0.0f, 0.61f, 0.56f, 0.64f, 0.0f, 0.0f, 0.45f,
0.48f, 0.53f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.73f, 0.61f, 0.63f, 0.0f, 0.0f, 0.85f, 0.66f, 0.78f, 0.0f, 0.0f, 0.80f, 0.64f,
0.83f, 0.0f, 0.0f, 0.61f, 0.59f, 0.71f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.60f, 0.40f, 0.53f, 0.0f, 0.0f, 0.63f, 0.29f, 0.53f,
0.0f, 0.0f, 0.57f, 0.25f, 0.55f, 0.0f, 0.0f, 0.48f, 0.32f, 0.56f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.72f, 0.60f, 0.61f, 0.0f,
0.0f, 0.84f, 0.64f, 0.76f, 0.0f, 0.0f, 0.78f, 0.62f, 0.81f, 0.0f, 0.0f, 0.60f, 0.57f,
0.70f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.52f, 0.46f, 0.44f, 0.0f, 0.0f, 0.63f, 0.54f, 0.59f, 0.0f, 0.0f, 0.59f, 0.54f, 0.63f,
0.0f, 0.0f, 0.44f, 0.46f, 0.52f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f
};
dataSet.AddPointField("pointvar", pointvar, nVerts);
return dataSet;
}
inline vtkm::cont::DataSet MakeTestDataSet::Make2DRectilinearDataSet0()
{
vtkm::cont::DataSetBuilderRectilinear dsb;

25
vtkm/filter/ContourTreeUniformAugmented.hxx
View File

@ -57,7 +57,7 @@
#include <vtkm/filter/ContourTreeUniformAugmented.h>
#include <vtkm/worklet/ContourTreeUniformAugmented.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
@ -161,14 +161,12 @@ vtkm::cont::DataSet ContourTreeAugmented::DoExecute(
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
}
// Use the GetRowsColsSlices struct defined in the header to collect the nRows, nCols, and nSlices information
// Use the GetPointDimensions struct defined in the header to collect the meshSize information
vtkm::worklet::ContourTreeAugmented worklet;
vtkm::Id nRows;
vtkm::Id nCols;
vtkm::Id nSlices = 1;
vtkm::Id3 meshSize;
const auto& cells = input.GetCellSet();
vtkm::filter::ApplyPolicyCellSet(cells, policy, *this)
.CastAndCall(vtkm::worklet::contourtree_augmented::GetRowsColsSlices(), nRows, nCols, nSlices);
.CastAndCall(vtkm::worklet::contourtree_augmented::GetPointDimensions(), meshSize);
// TODO blockIndex needs to change if we have multiple blocks per MPI rank and DoExecute is called for multiple blocks
std::size_t blockIndex = 0;
@ -193,9 +191,7 @@ vtkm::cont::DataSet ContourTreeAugmented::DoExecute(
MultiBlockTreeHelper ? MultiBlockTreeHelper->LocalSortOrders[blockIndex]
: this->MeshSortOrder,
this->NumIterations,
nRows,
nCols,
nSlices,
meshSize,
this->UseMarchingCubes,
compRegularStruct);
@ -319,7 +315,7 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(
// create the local data block structure
localDataBlocks[bi] = new vtkm::worklet::contourtree_distributed::ContourTreeBlockData<T>();
localDataBlocks[bi]->NumVertices = currContourTreeMesh->NumVertices;
localDataBlocks[bi]->SortOrder = currContourTreeMesh->SortOrder;
// localDataBlocks[bi]->SortOrder = currContourTreeMesh->SortOrder;
localDataBlocks[bi]->SortedValue = currContourTreeMesh->SortedValues;
localDataBlocks[bi]->GlobalMeshIndex = currContourTreeMesh->GlobalMeshIndex;
localDataBlocks[bi]->Neighbours = currContourTreeMesh->Neighbours;
@ -411,7 +407,9 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(
// Construct the contour tree mesh from the last block
vtkm::worklet::contourtree_augmented::ContourTreeMesh<T> contourTreeMeshOut;
contourTreeMeshOut.NumVertices = localDataBlocks[0]->NumVertices;
contourTreeMeshOut.SortOrder = localDataBlocks[0]->SortOrder;
// contourTreeMeshOut.SortOrder = localDataBlocks[0]->SortOrder;
contourTreeMeshOut.SortOrder = vtkm::cont::ArrayHandleIndex(contourTreeMeshOut.NumVertices);
contourTreeMeshOut.SortIndices = vtkm::cont::ArrayHandleIndex(contourTreeMeshOut.NumVertices);
contourTreeMeshOut.SortedValues = localDataBlocks[0]->SortedValue;
contourTreeMeshOut.GlobalMeshIndex = localDataBlocks[0]->GlobalMeshIndex;
contourTreeMeshOut.Neighbours = localDataBlocks[0]->Neighbours;
@ -424,10 +422,7 @@ VTKM_CONT void ContourTreeAugmented::DoPostExecute(
maxIdx[1] = maxIdx[1] - 1;
maxIdx[2] = maxIdx[2] > 0 ? (maxIdx[2] - 1) : 0;
auto meshBoundaryExecObj = contourTreeMeshOut.GetMeshBoundaryExecutionObject(
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[0],
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize[1],
minIdx,
maxIdx);
this->MultiBlockTreeHelper->MultiBlockSpatialDecomposition.GlobalSize, minIdx, maxIdx);
// Run the worklet to compute the final contour tree
worklet.Run(
contourTreeMeshOut.SortedValues, // Unused param. Provide something to keep API happy

245
vtkm/filter/testing/UnitTestContourTreeUniformAugmentedFilter.cxx
View File

@ -78,19 +78,29 @@ private:
//
// Internal helper function to execute the contour tree and save repeat code in tests
//
// datSets: 0 -> 5x5.txt (2D), 1 -> 8x9test.txt (2D), 2-> 5b.txt (3D)
vtkm::filter::ContourTreeAugmented RunContourTree(bool useMarchingCubes,
unsigned int computeRegularStructure,
unsigned int numDims = 2) const
unsigned int dataSetNo) const
{
// Create the input uniform cell set with values to contour
vtkm::cont::DataSet dataSet;
if (numDims == 2)
switch (dataSetNo)
{
dataSet = MakeTestDataSet().Make2DUniformDataSet1();
}
else if (numDims == 3)
{
dataSet = MakeTestDataSet().Make3DUniformDataSet1();
case 0:
dataSet = MakeTestDataSet().Make2DUniformDataSet1();
break;
case 1:
dataSet = MakeTestDataSet().Make2DUniformDataSet3();
break;
case 2:
dataSet = MakeTestDataSet().Make3DUniformDataSet1();
break;
case 3:
dataSet = MakeTestDataSet().Make3DUniformDataSet4();
break;
default:
VTKM_TEST_ASSERT(false);
}
vtkm::filter::ContourTreeAugmented filter(useMarchingCubes, computeRegularStructure);
filter.SetActiveField("pointvar");
@ -102,14 +112,15 @@ public:
//
// Create a uniform 2D structured cell set as input with values for contours
//
void TestContourTree_Mesh2D_Freudenthal(unsigned int computeRegularStructure = 1) const
void TestContourTree_Mesh2D_Freudenthal_SquareExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 2D Mesh. computeRegularStructure="
<< computeRegularStructure << std::endl;
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(false, // no marching cubes,
computeRegularStructure, // compute regular structure
2 // use 2D mesh
0 // use 5x5.txt
);
// Compute the saddle peaks to make sure the contour tree is correct
@ -148,7 +159,58 @@ public:
"Wrong result for ContourTree filter");
}
void TestContourTree_Mesh3D_Freudenthal(unsigned int computeRegularStructure = 1) const
void TestContourTree_Mesh2D_Freudenthal_NonSquareExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 2D Mesh. computeRegularStructure="
<< computeRegularStructure << std::endl;
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(false, // no marching cubes,
computeRegularStructure, // compute regular structure
1 // use 8x9test.txt
);
// Compute the saddle peaks to make sure the contour tree is correct
vtkm::worklet::contourtree_augmented::EdgePairArray saddlePeak;
vtkm::worklet::contourtree_augmented::ProcessContourTree::CollectSortedSuperarcs(
filter.GetContourTree(), filter.GetSortOrder(), saddlePeak);
// Print the contour tree we computed
std::cout << "Computed Contour Tree" << std::endl;
vtkm::worklet::contourtree_augmented::PrintEdgePairArray(saddlePeak);
// Print the expected contour tree
std::cout << "Expected Contour Tree" << std::endl;
std::cout << " 10 20" << std::endl;
std::cout << " 20 34" << std::endl;
std::cout << " 20 38" << std::endl;
std::cout << " 20 61" << std::endl;
std::cout << " 23 34" << std::endl;
std::cout << " 24 34" << std::endl;
std::cout << " 50 61" << std::endl;
std::cout << " 61 71" << std::endl;
VTKM_TEST_ASSERT(test_equal(saddlePeak.GetNumberOfValues(), 8),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(0), vtkm::make_Pair(10, 20)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(1), vtkm::make_Pair(20, 34)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(2), vtkm::make_Pair(20, 38)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(3), vtkm::make_Pair(20, 61)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(4), vtkm::make_Pair(23, 34)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(5), vtkm::make_Pair(24, 34)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(6), vtkm::make_Pair(50, 61)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(7), vtkm::make_Pair(61, 71)),
"Wrong result for ContourTree filter");
}
void TestContourTree_Mesh3D_Freudenthal_CubicExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 3D Mesh. computeRegularStructure="
<< computeRegularStructure << std::endl;
@ -157,7 +219,7 @@ public:
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(false, // no marching cubes,
computeRegularStructure, // compute regular structure
3 // use 2D mesh
2 // use 5b.txt (3D) mesh
);
// Compute the saddle peaks to make sure the contour tree is correct
@ -203,7 +265,64 @@ public:
"Wrong result for ContourTree filter");
}
void TestContourTree_Mesh3D_MarchingCubes(unsigned int computeRegularStructure = 1) const
void TestContourTree_Mesh3D_Freudenthal_NonCubicExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 3D Mesh. computeRegularStructure="
<< computeRegularStructure << std::endl;
// Execute the filter
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(false, // no marching cubes,
computeRegularStructure, // compute regular structure
3 // use 5b.txt (3D) upsampled to 5x6x7 mesh
);
// Compute the saddle peaks to make sure the contour tree is correct
vtkm::worklet::contourtree_augmented::EdgePairArray saddlePeak;
vtkm::worklet::contourtree_augmented::ProcessContourTree::CollectSortedSuperarcs(
filter.GetContourTree(), filter.GetSortOrder(), saddlePeak);
// Print the contour tree we computed
std::cout << "Computed Contour Tree" << std::endl;
vtkm::worklet::contourtree_augmented::PrintEdgePairArray(saddlePeak);
// Print the expected contour tree
std::cout << "Expected Contour Tree" << std::endl;
std::cout << " 0 112" << std::endl;
std::cout << " 71 72" << std::endl;
std::cout << " 72 78" << std::endl;
std::cout << " 72 101" << std::endl;
std::cout << " 101 112" << std::endl;
std::cout << " 101 132" << std::endl;
std::cout << " 107 112" << std::endl;
std::cout << " 131 132" << std::endl;
std::cout << " 132 138" << std::endl;
// Make sure the contour tree is correct
VTKM_TEST_ASSERT(test_equal(saddlePeak.GetNumberOfValues(), 9),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(0), vtkm::make_Pair(0, 112)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(1), vtkm::make_Pair(71, 72)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(2), vtkm::make_Pair(72, 78)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(3), vtkm::make_Pair(72, 101)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(4), vtkm::make_Pair(101, 112)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(5), vtkm::make_Pair(101, 132)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(6), vtkm::make_Pair(107, 112)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(7), vtkm::make_Pair(131, 132)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(8), vtkm::make_Pair(132, 138)),
"Wrong result for ContourTree filter");
}
void TestContourTree_Mesh3D_MarchingCubes_CubicExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 3D Mesh Marching Cubes. computeRegularStructure="
<< computeRegularStructure << std::endl;
@ -212,7 +331,7 @@ public:
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(true, // no marching cubes,
computeRegularStructure, // compute regular structure
3 // use 3D mesh
2 // use 5b.txt (3D) mesh
);
// Compute the saddle peaks to make sure the contour tree is correct
@ -263,28 +382,110 @@ public:
"Wrong result for ContourTree filter");
}
void TestContourTree_Mesh3D_MarchingCubes_NonCubicExtents(
unsigned int computeRegularStructure = 1) const
{
std::cout << "Testing ContourTree_Augmented 3D Mesh Marching Cubes. computeRegularStructure="
<< computeRegularStructure << std::endl;
// Execute the filter
vtkm::filter::ContourTreeAugmented filter =
RunContourTree(true, // no marching cubes,
computeRegularStructure, // compute regular structure
3 // use 5b.txt (3D) upsampled to 5x6x7 mesh
);
// Compute the saddle peaks to make sure the contour tree is correct
vtkm::worklet::contourtree_augmented::EdgePairArray saddlePeak;
vtkm::worklet::contourtree_augmented::ProcessContourTree::CollectSortedSuperarcs(
filter.GetContourTree(), filter.GetSortOrder(), saddlePeak);
// Print the contour tree we computed
std::cout << "Computed Contour Tree" << std::endl;
vtkm::worklet::contourtree_augmented::PrintEdgePairArray(saddlePeak);
// Print the expected contour tree
std::cout << "Expected Contour Tree" << std::endl;
std::cout << " 0 203" << std::endl;
std::cout << " 71 72" << std::endl;
std::cout << " 72 78" << std::endl;
std::cout << " 72 101" << std::endl;
std::cout << " 101 112" << std::endl;
std::cout << " 101 132" << std::endl;
std::cout << " 107 112" << std::endl;
std::cout << " 112 203" << std::endl;
std::cout << " 131 132" << std::endl;
std::cout << " 132 138" << std::endl;
std::cout << " 203 209" << std::endl;
VTKM_TEST_ASSERT(test_equal(saddlePeak.GetNumberOfValues(), 11),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(0), vtkm::make_Pair(0, 203)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(1), vtkm::make_Pair(71, 72)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(2), vtkm::make_Pair(72, 78)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(3), vtkm::make_Pair(72, 101)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(4), vtkm::make_Pair(101, 112)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(5), vtkm::make_Pair(101, 132)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(6), vtkm::make_Pair(107, 112)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(7), vtkm::make_Pair(112, 203)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(8), vtkm::make_Pair(131, 132)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(9), vtkm::make_Pair(132, 138)),
"Wrong result for ContourTree filter");
VTKM_TEST_ASSERT(test_equal(saddlePeak.WritePortal().Get(10), vtkm::make_Pair(203, 209)),
"Wrong result for ContourTree filter");
}
void operator()() const
{
// Test 2D Freudenthal with augmentation
this->TestContourTree_Mesh2D_Freudenthal(1);
this->TestContourTree_Mesh2D_Freudenthal_SquareExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh2D_Freudenthal(0);
this->TestContourTree_Mesh2D_Freudenthal_SquareExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh2D_Freudenthal(2);
this->TestContourTree_Mesh2D_Freudenthal_SquareExtents(2);
// Test 2D Freudenthal with augmentation
this->TestContourTree_Mesh2D_Freudenthal_NonSquareExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh2D_Freudenthal_NonSquareExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh2D_Freudenthal_NonSquareExtents(2);
// Test 3D Freudenthal with augmentation
this->TestContourTree_Mesh3D_Freudenthal(1);
this->TestContourTree_Mesh3D_Freudenthal_CubicExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh3D_Freudenthal(0);
this->TestContourTree_Mesh3D_Freudenthal_CubicExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh3D_Freudenthal(2);
this->TestContourTree_Mesh3D_Freudenthal_CubicExtents(2);
// Test 3D Freudenthal with augmentation
this->TestContourTree_Mesh3D_Freudenthal_NonCubicExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh3D_Freudenthal_NonCubicExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh3D_Freudenthal_NonCubicExtents(2);
// Test 3D marching cubes with augmentation
this->TestContourTree_Mesh3D_MarchingCubes(1);
this->TestContourTree_Mesh3D_MarchingCubes_CubicExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh3D_MarchingCubes(0);
this->TestContourTree_Mesh3D_MarchingCubes_CubicExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh3D_MarchingCubes(2);
this->TestContourTree_Mesh3D_MarchingCubes_CubicExtents(2);
// Test 3D marching cubes with augmentation
this->TestContourTree_Mesh3D_MarchingCubes_NonCubicExtents(1);
// Make sure the contour tree does not change when we disable augmentation
this->TestContourTree_Mesh3D_MarchingCubes_NonCubicExtents(0);
// Make sure the contour tree does not change when we use boundary augmentation
this->TestContourTree_Mesh3D_MarchingCubes_NonCubicExtents(2);
}
};
}

41
vtkm/worklet/ContourTreeUniformAugmented.h
View File

@ -71,14 +71,14 @@
#include <vtkm/worklet/contourtree_augmented/ActiveGraph.h>
#include <vtkm/worklet/contourtree_augmented/ContourTree.h>
#include <vtkm/worklet/contourtree_augmented/ContourTreeMaker.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/MergeTree.h>
#include <vtkm/worklet/contourtree_augmented/MeshExtrema.h>
#include <vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary2D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary2D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
namespace vtkm
{
@ -109,15 +109,17 @@ public:
* the full regular strucuture this is not needed because all vertices
* (including the boundary) will be addded to the tree anyways.
*/
template <typename FieldType, typename StorageType>
void Run(const vtkm::cont::ArrayHandle<FieldType, StorageType>
fieldArray, // TODO: We really should not need this
contourtree_augmented::ContourTreeMesh<FieldType>& mesh,
template <typename FieldType,
typename StorageType,
typename MeshType,
typename MeshBoundaryMeshExecType>
void Run(const vtkm::cont::ArrayHandle<FieldType, StorageType> fieldArray,
MeshType& mesh,
contourtree_augmented::ContourTree& contourTree,
contourtree_augmented::IdArrayType sortOrder,
contourtree_augmented::IdArrayType& sortOrder,
vtkm::Id& nIterations,
unsigned int computeRegularStructure,
const contourtree_augmented::MeshBoundaryContourTreeMeshExec& meshBoundary)
const MeshBoundaryMeshExecType& meshBoundary)
{
RunContourTree(
fieldArray, // Just a place-holder to fill the required field. Used when calling SortData on the contour tree which is a no-op
@ -157,18 +159,16 @@ public:
contourtree_augmented::ContourTree& contourTree,
contourtree_augmented::IdArrayType& sortOrder,
vtkm::Id& nIterations,
const vtkm::Id nRows,
const vtkm::Id nCols,
const vtkm::Id nSlices = 1,
const vtkm::Id3 meshSize,
bool useMarchingCubes = false,
unsigned int computeRegularStructure = 1)
{
using namespace vtkm::worklet::contourtree_augmented;
// 2D Contour Tree
if (nSlices == 1)
if (meshSize[2] == 1)
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_2D_Freudenthal<FieldType, StorageType> mesh(nRows, nCols);
DataSetMeshTriangulation2DFreudenthal mesh(vtkm::Id2{ meshSize[0], meshSize[1] });
// Run the contour tree on the mesh
RunContourTree(fieldArray,
contourTree,
@ -183,7 +183,7 @@ public:
else if (useMarchingCubes)
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_3D_MarchingCubes<FieldType, StorageType> mesh(nRows, nCols, nSlices);
DataSetMeshTriangulation3DMarchingCubes mesh(meshSize);
// Run the contour tree on the mesh
RunContourTree(fieldArray,
contourTree,
@ -198,7 +198,7 @@ public:
else
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_3D_Freudenthal<FieldType, StorageType> mesh(nRows, nCols, nSlices);
DataSetMeshTriangulation3DFreudenthal mesh(meshSize);
// Run the contour tree on the mesh
RunContourTree(fieldArray,
contourTree,
@ -355,7 +355,10 @@ private:
// Collect the output data
nIterations = treeMaker.ContourTreeResult.NumIterations;
sortOrder = mesh.SortOrder;
// Need to make a copy of sortOrder since ContourTreeMesh uses a smart array handle
// TODO: Check if we can just make sortOrder a return array with variable type or if we can make the SortOrder return optional
vtkm::cont::Algorithm::Copy(mesh.SortOrder, sortOrder);
// sortOrder = mesh.SortOrder;
// ProcessContourTree::CollectSortedSuperarcs<DeviceAdapter>(contourTree, mesh.SortOrder, saddlePeak);
// contourTree.SortedArcPrint(mesh.SortOrder);
// contourTree.PrintDotSuperStructure();

6
vtkm/worklet/contourtree_augmented/CMakeLists.txt
View File

@ -13,8 +13,8 @@ set(headers
ActiveGraph.h
ContourTree.h
ContourTreeMaker.h
DataSetMesh.h
MergeTree.h
Mesh_DEM_Triangulation.h
MeshExtrema.h
PointerDoubling.h
PrintVectors.h
@ -22,8 +22,8 @@ set(headers
Types.h
)
#-----------------------------------------------------------------------------
add_subdirectory(mesh_dem)
add_subdirectory(mesh_dem_meshtypes)
add_subdirectory(data_set_mesh)
add_subdirectory(meshtypes)
add_subdirectory(meshextrema)
add_subdirectory(activegraph)
add_subdirectory(contourtreemaker)

28
vtkm/worklet/contourtree_augmented/ContourTreeMaker.h
View File

@ -58,9 +58,9 @@
// local includes
#include <vtkm/worklet/contourtree_augmented/ArrayTransforms.h>
#include <vtkm/worklet/contourtree_augmented/ContourTree.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/MergeTree.h>
#include <vtkm/worklet/contourtree_augmented/MeshExtrema.h>
#include <vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h>
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
@ -534,20 +534,26 @@ void ContourTreeMaker::ComputeBoundaryRegularStructure(
// for (indexType supernode = 0; supernode < this->ContourTreeResult.Supernodes.size(); supernode++)
// superparents[this->ContourTreeResult.Supernodes[supernode]] = supernode;
IdArrayType sortOrder; // TODO/FIX: Why is this copy necessary?
vtkm::cont::Algorithm::Copy(mesh.SortOrder, sortOrder);
// Second step - for all remaining (regular) nodes, locate the superarc to which they belong
contourtree_maker_inc_ns::ComputeRegularStructure_LocateSuperarcsOnBoundary
locateSuperarcsOnBoundaryWorklet(this->ContourTreeResult.Hypernodes.GetNumberOfValues(),
this->ContourTreeResult.Supernodes.GetNumberOfValues());
this->Invoke(locateSuperarcsOnBoundaryWorklet,
superparents, // (input/output)
this->ContourTreeResult.WhenTransferred, // (input)
this->ContourTreeResult.Hyperparents, // (input)
this->ContourTreeResult.Hyperarcs, // (input)
this->ContourTreeResult.Hypernodes, // (input)
this->ContourTreeResult.Supernodes, // (input)
meshExtrema.Peaks, // (input)
meshExtrema.Pits, // (input)
meshBoundaryExecObj); // (input)
this->Invoke(
locateSuperarcsOnBoundaryWorklet,
superparents, // (input/output)
this->ContourTreeResult.WhenTransferred, // (input)
this->ContourTreeResult.Hyperparents, // (input)
this->ContourTreeResult.Hyperarcs, // (input)
this->ContourTreeResult.Hypernodes, // (input)
this->ContourTreeResult.Supernodes, // (input)
meshExtrema.Peaks, // (input)
meshExtrema.Pits, // (input)
sortOrder, // (input)
//mesh.SortOrder, // (input) // TODO/FIXME: Why can't I just pass this?
meshBoundaryExecObj); // (input)
// We have now set the superparent correctly for each node, and need to sort them to get the correct regular arcs
// DAVID "ContourTreeMaker.h" line 338

292
vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h → vtkm/worklet/contourtree_augmented/DataSetMesh.h
View File

@ -68,10 +68,8 @@
//
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_h
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
@ -82,12 +80,11 @@
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/SimulatedSimplicityComperator.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/SortIndices.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/SimulatedSimplicityComperator.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/SortIndices.h>
//Define namespace alias for the freudenthal types to make the code a bit more readable
namespace mesh_dem_ns = vtkm::worklet::contourtree_augmented::mesh_dem;
namespace vtkm
{
@ -95,17 +92,15 @@ namespace worklet
{
namespace contourtree_augmented
{
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation
class DataSetMesh
{
public:
// common mesh size parameters
// common mesh size parameter, use all three dimensions ofr MeshSize with third determining if 2D or 3D
// (convention: MeshSize[2] is always >= 1, even for empty data set, so that we can detect 2D
// data as MeshSize[2] == 1)
vtkm::Id3 MeshSize;
vtkm::Id NumVertices, NumLogSteps;
// Define dimensionality of the mesh
vtkm::Id NumDims;
// Array with the sorted order of the mesh vertices
IdArrayType SortOrder;
@ -114,152 +109,146 @@ public:
IdArrayType SortIndices;
//empty constructor
Mesh_DEM_Triangulation()
: NumVertices(0)
, NumLogSteps(0)
, NumDims(2)
DataSetMesh()
: MeshSize{ 0, 0, 1 } // Always set third dimension to 1 for easy detection of 2D vs 3D
, NumVertices(0)
, NumLogSteps(1)
{
}
// base constructor
DataSetMesh(vtkm::Id3 meshSize)
: MeshSize{ meshSize }
, NumVertices{ meshSize[0] * meshSize[1] * meshSize[2] }
// per convention meshSize[2] == 1 for 2D
, NumLogSteps(1)
{
// Per convention the third dimension should be 1 (even for an empty
// mesh) or higher to make it easier to check for 2D vs. 3D data)
VTKM_ASSERT(MeshSize[2] >= 1);
// TODO/FIXME: An empty mesh will likely cause a crash down the
// road anyway, so we may want to detect that case and handle
// it appropriately.
// Compute the number of log-jumping steps (i.e. lg_2 (NumVertices))
// this->NumLogSteps = 1; // already set in initializer list
for (vtkm::Id shifter = this->NumVertices; shifter > 0; shifter >>= 1)
this->NumLogSteps++;
}
virtual ~DataSetMesh() {}
// Getter function for NumVertices
vtkm::Id GetNumberOfVertices() const { return this->NumVertices; }
// sorts the data and initializes the sortIndex & indexReverse
// Sorts the data and initializes SortOrder & SortIndices
template <typename T, typename StorageType>
void SortData(const vtkm::cont::ArrayHandle<T, StorageType>& values);
/// Routine to return the global IDs for a set of vertices
/// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
/// @param[in] meshIds Array with sort Ids to be converted from local to global Ids
/// @param[in] localToGlobalIdRelabeler This parameter is the IdRelabeler
/// used to transform local to global Ids. The relabeler relies on the
/// decomposition of the global mesh which is not know by this block.
inline vtkm::cont::ArrayHandleTransform<
vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType>,
mesh_dem::IdRelabeler>
GetGlobalIdsFromSortIndices(const IdArrayType& sortIds,
const mesh_dem::IdRelabeler* localToGlobalIdRelabeler) const
{ // GetGlobalIDsFromSortIndices()
auto permutedSortOrder = vtkm::cont::make_ArrayHandlePermutation(sortIds, this->SortOrder);
return vtkm::cont::ArrayHandleTransform<
vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType>,
mesh_dem::IdRelabeler>(permutedSortOrder, *localToGlobalIdRelabeler);
} // GetGlobalIDsFromSortIndices()
/// Routine to return the global IDs for a set of vertices
/// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
/// @param[in] meshIds Array with mesh Ids to be converted from local to global Ids
/// @param[in] localToGlobalIdRelabeler This parameter is the IdRelabeler
/// used to transform local to global Ids. The relabeler relies on the
/// decomposition of the global mesh which is not know by this block.
inline vtkm::cont::ArrayHandleTransform<IdArrayType, mesh_dem::IdRelabeler>
GetGlobalIdsFromMeshIndices(const IdArrayType& meshIds,
const mesh_dem::IdRelabeler* localToGlobalIdRelabeler) const
{ // GetGlobalIDsFromMeshIndices()
return vtkm::cont::ArrayHandleTransform<IdArrayType, mesh_dem::IdRelabeler>(
meshIds, *localToGlobalIdRelabeler);
} // GetGlobalIDsFromMeshIndices()
//routine that dumps out the contents of the mesh
void DebugPrint(const char* message, const char* fileName, long lineNum);
protected:
virtual void DebugPrintExtends() = 0;
virtual void DebugPrintValues(const vtkm::cont::ArrayHandle<T, StorageType>& values) = 0;
}; // class Mesh_DEM_Triangulation
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation_2D : public Mesh_DEM_Triangulation<T, StorageType>
{
public:
// 2D mesh size parameters
vtkm::Id NumColumns, NumRows;
// Maximum outdegree
static constexpr int MAX_OUTDEGREE = 3;
// empty constructor
Mesh_DEM_Triangulation_2D()
: Mesh_DEM_Triangulation<T, StorageType>()
, NumColumns(0)
, NumRows(0)
{
this->NumDims = 2;
}
// base constructor
Mesh_DEM_Triangulation_2D(vtkm::Id ncols, vtkm::Id nrows)
: Mesh_DEM_Triangulation<T, StorageType>()
, NumColumns(ncols)
, NumRows(nrows)
{
this->NumDims = 2;
this->NumVertices = NumRows * NumColumns;
// compute the number of log-jumping steps (i.e. lg_2 (NumVertices))
this->NumLogSteps = 1;
for (vtkm::Id shifter = this->NumVertices; shifter > 0; shifter >>= 1)
this->NumLogSteps++;
}
protected:
virtual void DebugPrintExtends();
virtual void DebugPrintValues(const vtkm::cont::ArrayHandle<T, StorageType>& values);
}; // class Mesh_DEM_Triangulation_2D
template <typename T, typename StorageType>
void DebugPrintValues(const vtkm::cont::ArrayHandle<T, StorageType>& values);
}; // class DataSetMesh
// Sorts the data and initialises the SortIndices & SortOrder
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation_3D : public Mesh_DEM_Triangulation<T, StorageType>
{
public:
// 2D mesh size parameters
vtkm::Id NumColumns, NumRows, NumSlices;
// Maximum outdegree
static constexpr int MAX_OUTDEGREE = 6; // True for Freudenthal and Marching Cubes
// empty constructor
Mesh_DEM_Triangulation_3D()
: Mesh_DEM_Triangulation<T, StorageType>()
, NumColumns(0)
, NumRows(0)
, NumSlices(0)
{
this->NumDims = 3;
}
// base constructor
Mesh_DEM_Triangulation_3D(vtkm::Id ncols, vtkm::Id nrows, vtkm::Id nslices)
: Mesh_DEM_Triangulation<T, StorageType>()
, NumColumns(ncols)
, NumRows(nrows)
, NumSlices(nslices)
{
this->NumDims = 3;
this->NumVertices = NumRows * NumColumns * NumSlices;
// compute the number of log-jumping steps (i.e. lg_2 (NumVertices))
this->NumLogSteps = 1;
for (vtkm::Id shifter = this->NumVertices; shifter > 0; shifter >>= 1)
this->NumLogSteps++;
}
protected:
virtual void DebugPrintExtends();
virtual void DebugPrintValues(const vtkm::cont::ArrayHandle<T, StorageType>& values);
}; // class Mesh_DEM_Triangulation_3D
// sorts the data and initialises the SortIndices & SortOrder
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation<T, StorageType>::SortData(
const vtkm::cont::ArrayHandle<T, StorageType>& values)
void DataSetMesh::SortData(const vtkm::cont::ArrayHandle<T, StorageType>& values)
{
// Define namespace alias for mesh dem worklets
namespace mesh_dem_worklets = vtkm::worklet::contourtree_augmented::mesh_dem;
// Make sure that the values have the correct size
assert(values.GetNumberOfValues() == NumVertices);
VTKM_ASSERT(values.GetNumberOfValues() == this->NumVertices);
// Make sure that we are not running on an empty mesh
VTKM_ASSERT(this->NumVertices > 0);
// Just in case, make sure that everything is cleaned up
SortIndices.ReleaseResources();
SortOrder.ReleaseResources();
this->SortIndices.ReleaseResources();
this->SortOrder.ReleaseResources();
// allocate memory for the sort arrays
SortOrder.Allocate(NumVertices);
SortIndices.Allocate(NumVertices);
SortOrder.Allocate(this->NumVertices);
SortIndices.Allocate(this->NumVertices);
// now sort the sort order vector by the values, i.e,. initialize the SortOrder member variable
vtkm::cont::ArrayHandleIndex initVertexIds(NumVertices); // create sequence 0, 1, .. NumVertices
vtkm::cont::ArrayCopy(initVertexIds, SortOrder);
vtkm::cont::ArrayHandleIndex initVertexIds(
this->NumVertices); // create sequence 0, 1, .. NumVertices
vtkm::cont::ArrayCopy(initVertexIds, this->SortOrder);
vtkm::cont::Algorithm::Sort(SortOrder,
vtkm::cont::Algorithm::Sort(this->SortOrder,
mesh_dem::SimulatedSimplicityIndexComparator<T, StorageType>(values));
// now set the index lookup, i.e., initialize the SortIndices member variable
// In serial this would be
// for (indexType vertex = 0; vertex < NumVertices; vertex++)
// SortIndices[SortOrder[vertex]] = vertex;
mesh_dem_worklets::SortIndices sortIndicesWorklet;
data_set_mesh::SortIndices sortIndicesWorklet;
vtkm::cont::Invoker invoke;
invoke(sortIndicesWorklet, SortOrder, SortIndices);
invoke(sortIndicesWorklet, this->SortOrder, this->SortIndices);
// Debug print statement
DebugPrint("Data Sorted", __FILE__, __LINE__);
DebugPrintValues(values);
} // SortData()
// Print mesh extends
void DataSetMesh::DebugPrintExtends()
{
// For compatibility with the output of the original PPP Implementation, print size
// as NumRows, NumColumns and NumSlices (if applicable)
PrintLabel("NumRows");
PrintIndexType(this->MeshSize[1]);
std::cout << std::endl;
PrintLabel("NumColumns");
PrintIndexType(this->MeshSize[0]);
std::cout << std::endl;
if (MeshSize[2] > 1)
{
PrintLabel("NumSlices");
PrintIndexType(this->MeshSize[2]);
std::cout << std::endl;
}
} // DebugPrintExtends
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation<T, StorageType>::DebugPrint(const char* message,
const char* fileName,
long lineNum)
void DataSetMesh::DebugPrint(const char* message, const char* fileName, long lineNum)
{ // DebugPrint()
#ifdef DEBUG_PRINT
std::cout << "------------------------------------------------------" << std::endl;
@ -270,13 +259,13 @@ void Mesh_DEM_Triangulation<T, StorageType>::DebugPrint(const char* message,
std::cout << "------------------------------------------------------" << std::endl;
//DebugPrintExtents();
PrintLabel("NumVertices");
PrintIndexType(NumVertices);
PrintIndexType(this->NumVertices);
std::cout << std::endl;
PrintLabel("NumLogSteps");
PrintIndexType(this->NumLogSteps);
std::cout << std::endl;
PrintIndices("Sort Indices", SortIndices);
PrintIndices("Sort Order", SortOrder);
PrintIndices("Sort Indices", this->SortIndices);
PrintIndices("Sort Order", this->SortOrder);
std::cout << std::endl;
#else
// Avoid unused parameter warning
@ -286,41 +275,13 @@ void Mesh_DEM_Triangulation<T, StorageType>::DebugPrint(const char* message,
#endif
} // DebugPrint()
// print mesh extends for 2D mesh
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_2D<T, StorageType>::DebugPrintExtends()
{
PrintLabel("NumRows");
PrintIndexType(NumRows);
std::cout << std::endl;
PrintLabel("NumColumns");
PrintIndexType(NumColumns);
std::cout << std::endl;
} // DebugPrintExtends for 2D
// print mesh extends for 3D mesh
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D<T, StorageType>::DebugPrintExtends()
{
PrintLabel("NumRows");
PrintIndexType(NumRows);
std::cout << std::endl;
PrintLabel("NumColumns");
PrintIndexType(NumColumns);
std::cout << std::endl;
PrintLabel("NumSlices");
PrintIndexType(NumSlices);
std::cout << std::endl;
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_2D<T, StorageType>::DebugPrintValues(
const vtkm::cont::ArrayHandle<T, StorageType>& values)
void DataSetMesh::DebugPrintValues(const vtkm::cont::ArrayHandle<T, StorageType>& values)
{
#ifdef DEBUG_PRINT
if (NumColumns > 0)
if (MeshSize[0] > 0)
{
PrintLabelledDataBlock<T, StorageType>("Value", values, NumColumns);
PrintLabelledDataBlock<T, StorageType>("Value", values, MeshSize[0]);
PrintSortedValues("Sorted Values", values, this->SortOrder);
}
PrintHeader(values.GetNumberOfValues());
@ -330,28 +291,13 @@ void Mesh_DEM_Triangulation_2D<T, StorageType>::DebugPrintValues(
#endif
} // DebugPrintValues
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D<T, StorageType>::DebugPrintValues(
const vtkm::cont::ArrayHandle<T, StorageType>& values)
{
#ifdef DEBUG_PRINT
if (NumColumns > 0)
{
PrintLabelledDataBlock<T, StorageType>("Value", values, NumColumns);
}
PrintHeader(values.GetNumberOfValues());
#else
// Avoid unused parameter warning
(void)values;
#endif
} // DebugPrintValues
} // namespace contourtree_augmented
} // worklet
} // vtkm
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/Freudenthal_2D_Triangulation.h> // include Mesh_DEM_Triangulation_2D_Freudenthal
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/Freudenthal_3D_Triangulation.h> // include Mesh_DEM_Triangulation_3D_Freudenthal
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MarchingCubes_3D_Triangulation.h> // include Mesh_DEM_Triangulation_3D_MarchinCubes
// Include specialized mesh classes providing triangulation/connectivity information
#include <vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation2DFreudenthal.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DFreudenthal.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DMarchingCubes.h>
#endif

15
vtkm/worklet/contourtree_augmented/MergeTree.h
View File

@ -59,7 +59,7 @@
// local includes
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
//VTKM includes
@ -250,16 +250,17 @@ inline void MergeTree::DebugPrintTree(const char* message,
{
std::cout << mesh.SortOrder.ReadPortal().Get(arc) << std::endl;
}
if (mesh.NumDims == 2)
{
if ((entry % mesh.NumColumns) == (mesh.NumColumns - 1))
if (mesh.MeshSize[2] == 1)
{ // 2D Mesh
if ((entry % mesh.MeshSize[0]) == (mesh.MeshSize[0] - 1))
{
std::cout << std::endl;
}
}
else if (mesh.NumDims == 3)
{
if ((entry % (mesh.NumColumns * mesh.NumRows)) == (mesh.NumColumns * mesh.NumRows - 1))
else
{ // 3D Mesh
if ((entry % (mesh.MeshSize[0] * mesh.MeshSize[1])) ==
(mesh.MeshSize[0] * mesh.MeshSize[1] - 1))
{
std::cout << std::endl;
}

2
vtkm/worklet/contourtree_augmented/MeshExtrema.h
View File

@ -65,7 +65,7 @@
#include <vtkm/worklet/contourtree_augmented/meshextrema/SetStarts.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/SortIndices.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/SortIndices.h>
namespace mesh_extrema_inc_ns = vtkm::worklet::contourtree_augmented::mesh_extrema_inc;

35
vtkm/worklet/contourtree_augmented/Types.h
View File

@ -249,44 +249,29 @@ public:
/// to determine the rows, cols, slices parameters from the
/// datasets so we can call the contour tree worklet properly.
///
struct GetRowsColsSlices
struct GetPointDimensions
{
//@{
/// Get the number of rows, cols, and slices of a vtkm::cont::CellSetStructured
/// @param[in] cells The input vtkm::cont::CellSetStructured
/// @param[out] nRows Number of rows (x) in the cell set
/// @param[out[ nCols Number of columns (y) in the cell set
/// @param[out] nSlices Number of slices (z) in the cell set
void operator()(const vtkm::cont::CellSetStructured<2>& cells,
vtkm::Id& nRows,
vtkm::Id& nCols,
vtkm::Id& nSlices) const
/// @param[out] pointDimensions mesh size (#cols, #rows #slices in old notation) with last dimension having a value of 1 for 2D data
void operator()(const vtkm::cont::CellSetStructured<2>& cells, vtkm::Id3& pointDimensions) const
{
vtkm::Id2 pointDimensions = cells.GetPointDimensions();
nRows = pointDimensions[0];
nCols = pointDimensions[1];
nSlices = 1;
vtkm::Id2 pointDimensions2D = cells.GetPointDimensions();
pointDimensions[0] = pointDimensions2D[0];
pointDimensions[1] = pointDimensions2D[1];
pointDimensions[2] = 1;
}
void operator()(const vtkm::cont::CellSetStructured<3>& cells,
vtkm::Id& nRows,
vtkm::Id& nCols,
vtkm::Id& nSlices) const
void operator()(const vtkm::cont::CellSetStructured<3>& cells, vtkm::Id3& pointDimensions) const
{
vtkm::Id3 pointDimensions = cells.GetPointDimensions();
nRows = pointDimensions[0];
nCols = pointDimensions[1];
nSlices = pointDimensions[2];
pointDimensions = cells.GetPointDimensions();
}
//@}
/// Raise ErrorBadValue if the input cell set is not a vtkm::cont::CellSetStructured<2> or <3>
template <typename T>
void operator()(const T& cells, vtkm::Id& nRows, vtkm::Id& nCols, vtkm::Id& nSlices) const
void operator()(const T&, vtkm::Id3&) const
{
(void)nRows;
(void)nCols;
(void)nSlices;
(void)cells;
throw vtkm::cont::ErrorBadValue("Expected 2D or 3D structured cell cet! ");
}
};

9
vtkm/worklet/contourtree_augmented/contourtreemaker/ComputeRegularStructure_LocateSuperarcs.h
View File

@ -415,16 +415,15 @@ public:
WholeArrayIn contourTreeSupernodes, // (input)
WholeArrayIn meshExtremaPeaks, // (input)
WholeArrayIn meshExtremaPits, // (input)
WholeArrayIn sortOrder, // (input)
ExecObject meshBoundary); // (input)
typedef void ExecutionSignature(_1, InputIndex, _2, _3, _4, _5, _6, _7, _8, _9);
typedef void ExecutionSignature(_1, InputIndex, _2, _3, _4, _5, _6, _7, _8, _9, _10);
using InputDomain = _1;
vtkm::Id NumHypernodes; // contourTree.Hypernodes.GetNumberOfValues()
vtkm::Id NumSupernodes; // contourTree.Supernodes.GetNumberOfValues()
vtkm::Id NumRows, NumColumns, NumSlices; // Mesh 2D or 3D - NumRows, NumColumns, NumSlices
// Default Constructor
VTKM_EXEC_CONT
ComputeRegularStructure_LocateSuperarcsOnBoundary(vtkm::Id numHypernodes, vtkm::Id numSupernodes)
@ -443,11 +442,13 @@ public:
const InFieldPortalType& contourTreeSupernodesPortal,
const InFieldPortalType& meshExtremaPeaksPortal,
const InFieldPortalType& meshExtremaPitsPortal,
const InFieldPortalType& sortOrderPortal,
const MeshBoundaryType& meshBoundary) const
{
// per node
// if the superparent is already set, it's a supernode, so skip it.
if (NoSuchElement(contourTreeSuperparentsPortal.Get(node)) && meshBoundary.liesOnBoundary(node))
if (NoSuchElement(contourTreeSuperparentsPortal.Get(node)) &&
meshBoundary.LiesOnBoundary(sortOrderPortal.Get(node)))
{ // regular nodes only
// we will need to prune top and bottom until one of them prunes past the node
vtkm::Id top = meshExtremaPeaksPortal.Get(node);

0
vtkm/worklet/contourtree_augmented/mesh_dem/CMakeLists.txt → vtkm/worklet/contourtree_augmented/data_set_mesh/CMakeLists.txt
View File

48
vtkm/worklet/contourtree_augmented/mesh_dem/IdRelabeler.h → vtkm/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h
View File

@ -75,57 +75,45 @@ namespace contourtree_augmented
namespace mesh_dem
{
/// A utility class that converts Ids from local to global given a mesh
class IdRelabeler
{
public:
VTKM_EXEC_CONT
IdRelabeler()
: InputStartRow(0)
, InputStartCol(0)
, InputStartSlice(0)
, InputNumRows(1)
, InputNumCols(1)
, OutputNumRows(1)
, OutputNumCol(1)
: LocalBlockOrigin{ 0, 0, 0 }
, LocalBlockSize{ 1, 1, 1 }
, GlobalSize{ 1, 1, 1 }
{
}
VTKM_EXEC_CONT
IdRelabeler(vtkm::Id iSR,
vtkm::Id iSC,
vtkm::Id iSS,
vtkm::Id iNR,
vtkm::Id iNC,
vtkm::Id oNR,
vtkm::Id oNC)
: InputStartRow(iSR)
, InputStartCol(iSC)
, InputStartSlice(iSS)
, InputNumRows(iNR)
, InputNumCols(iNC)
, OutputNumRows(oNR)
, OutputNumCol(oNC)
IdRelabeler(vtkm::Id3 lBO, vtkm::Id3 lBS, vtkm::Id3 gS)
: LocalBlockOrigin(lBO)
, LocalBlockSize(lBS)
, GlobalSize(gS)
{
}
VTKM_EXEC_CONT
vtkm::Id operator()(vtkm::Id v) const
{
vtkm::Id r = InputStartRow + ((v % (InputNumRows * InputNumCols)) / InputNumCols);
vtkm::Id c = InputStartCol + (v % InputNumCols);
vtkm::Id s = InputStartSlice + v / (InputNumRows * InputNumCols);
// Translate v into mesh coordinates and add offset
vtkm::Id3 pos{ LocalBlockOrigin[0] + (v % LocalBlockSize[0]),
LocalBlockOrigin[1] +
(v % (LocalBlockSize[1] * LocalBlockSize[0]) / LocalBlockSize[0]),
LocalBlockOrigin[2] + (v / (LocalBlockSize[1] * LocalBlockSize[0])) };
return (s * OutputNumRows + r) * OutputNumCol + c;
// Translate mesh coordinates into global Id
return (pos[2] * GlobalSize[1] + pos[1]) * GlobalSize[0] + pos[0];
}
private:
vtkm::Id InputStartRow, InputStartCol, InputStartSlice;
vtkm::Id InputNumRows, InputNumCols;
vtkm::Id OutputNumRows, OutputNumCol;
vtkm::Id3 LocalBlockOrigin;
vtkm::Id3 LocalBlockSize;
vtkm::Id3 GlobalSize;
};
} // namespace mesh_dem
} // namespace contourtree_augmented
} // namespace worklet

31
vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure2D.h → vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure2D.h
View File

@ -50,19 +50,18 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_execution_object_mesh_2d_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_execution_object_mesh_2d_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_execution_object_mesh_2d_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_execution_object_mesh_2d_h
#include <vtkm/Types.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace mesh_dem
namespace data_set_mesh
{
// Worklet for computing the sort indices from the sort order
@ -72,35 +71,31 @@ class MeshStructure2D
public:
VTKM_EXEC_CONT
MeshStructure2D()
: NumColumns(0)
, NumRows(0)
: MeshSize{ 0, 0 }
{
}
VTKM_EXEC_CONT
MeshStructure2D(vtkm::Id ncols, vtkm::Id nrows)
: NumColumns(ncols)
, NumRows(nrows)
MeshStructure2D(vtkm::Id2 meshSize)
: MeshSize(meshSize)
{
}
// number of mesh vertices
VTKM_EXEC_CONT
vtkm::Id GetNumberOfVertices() const { return (this->NumRows * this->NumColumns); }
vtkm::Id GetNumberOfVertices() const { return (this->MeshSize[0] * this->MeshSize[1]); }
// vertex row - integer divide by columns
VTKM_EXEC
inline vtkm::Id VertexRow(vtkm::Id v) const { return v / NumColumns; }
// verteck column -- integer modulus by columns
VTKM_EXEC
inline vtkm::Id VertexColumn(vtkm::Id v) const { return v % NumColumns; }
inline vtkm::Id2 VertexPos(vtkm::Id v) const
{
return vtkm::Id2{ v % this->MeshSize[0], v / this->MeshSize[0] };
}
//vertex ID - row * ncols + col
VTKM_EXEC
inline vtkm::Id VertexId(vtkm::Id r, vtkm::Id c) const { return r * NumColumns + c; }
inline vtkm::Id VertexId(vtkm::Id2 pos) const { return pos[1] * this->MeshSize[0] + pos[0]; }
vtkm::Id NumColumns, NumRows;
vtkm::Id2 MeshSize;
}; // MeshStructure2D

43
vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure3D.h → vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure3D.h
View File

@ -50,19 +50,18 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_execution_object_mesh_3d_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_execution_object_mesh_3d_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_execution_object_mesh_3d_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_execution_object_mesh_3d_h
#include <vtkm/Types.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace mesh_dem
namespace data_set_mesh
{
// Worklet for computing the sort indices from the sort order
@ -72,17 +71,13 @@ class MeshStructure3D
public:
VTKM_EXEC_CONT
MeshStructure3D()
: NumColumns(0)
, NumRows(0)
, NumSlices(0)
: MeshSize{ 0, 0, 0 }
{
}
VTKM_EXEC_CONT
MeshStructure3D(vtkm::Id ncols, vtkm::Id nrows, vtkm::Id nslices)
: NumColumns(ncols)
, NumRows(nrows)
, NumSlices(nslices)
MeshStructure3D(vtkm::Id3 meshSize)
: MeshSize(meshSize)
{
}
@ -90,29 +85,27 @@ public:
VTKM_EXEC_CONT
vtkm::Id GetNumberOfVertices() const
{
return (this->NumRows * this->NumColumns * this->NumSlices);
return (this->MeshSize[0] * this->MeshSize[1] * this->MeshSize[2]);
}
// vertex row - integer modulus by (NumRows&NumColumns) and integer divide by columns
VTKM_EXEC
vtkm::Id VertexRow(vtkm::Id v) const { return (v % (NumRows * NumColumns)) / NumColumns; }
// vertex column -- integer modulus by columns
VTKM_EXEC
vtkm::Id VertexColumn(vtkm::Id v) const { return v % NumColumns; }
// vertex slice -- integer divide by (NumRows*NumColumns)
VTKM_EXEC
vtkm::Id VertexSlice(vtkm::Id v) const { return v / (NumRows * NumColumns); }
vtkm::Id3 VertexPos(vtkm::Id v) const
{
return vtkm::Id3{
v % this->MeshSize[0], // column
(v % (this->MeshSize[1] * this->MeshSize[0])) / this->MeshSize[0], // row
v / (this->MeshSize[1] * this->MeshSize[0]) // slice
};
}
//vertex ID - row * ncols + col
VTKM_EXEC
vtkm::Id VertexId(vtkm::Id s, vtkm::Id r, vtkm::Id c) const
vtkm::Id VertexId(vtkm::Id3 pos) const
{
return (s * NumRows + r) * NumColumns + c;
return (pos[2] * this->MeshSize[1] + pos[1]) * this->MeshSize[0] + pos[0];
}
vtkm::Id NumColumns, NumRows, NumSlices;
vtkm::Id3 MeshSize;
}; // Mesh_DEM_2D_ExecutionObject

0
vtkm/worklet/contourtree_augmented/mesh_dem/SimulatedSimplicityComperator.h → vtkm/worklet/contourtree_augmented/data_set_mesh/SimulatedSimplicityComperator.h
View File

4
vtkm/worklet/contourtree_augmented/mesh_dem/SortIndices.h → vtkm/worklet/contourtree_augmented/data_set_mesh/SortIndices.h
View File

@ -61,7 +61,7 @@ namespace worklet
{
namespace contourtree_augmented
{
namespace mesh_dem
namespace data_set_mesh
{
// Worklet for computing the sort indices from the sort order
@ -86,7 +86,7 @@ public:
}
}; // Mesh2D_DEM_VertexStarter
} // namespace mesh_dem_triangulation_worklets
} // namespace data_set_mesh
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm

154
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary3D.h
View File

@ -1,154 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt 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) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
// This header contains a collection of classes used to describe the boundary
// of a mesh, for each main mesh type (i.e., 2D, 3D, and ContourTreeMesh).
// For each mesh type, there are two classes, the actual boundary desriptor
// class and an ExectionObject class with the PrepareForInput function that
// VTKm expects to generate the object for the execution environment.
#ifndef vtk_m_worklet_contourtree_augmented_mesh_boundary_mesh_boundary_3d_h
#define vtk_m_worklet_contourtree_augmented_mesh_boundary_mesh_boundary_3d_h
#include <cstdlib>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure3D.h>
#include <vtkm/cont/ExecutionObjectBase.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
template <typename DeviceTag>
class MeshBoundary3D : public vtkm::cont::ExecutionObjectBase
{
public:
// Sort indicies types
using SortOrderPortalType = typename IdArrayType::template ExecutionTypes<DeviceTag>::PortalConst;
VTKM_EXEC_CONT
MeshBoundary3D()
: MeshStructure(mesh_dem::MeshStructure3D<DeviceTag>(0, 0, 0))
{
}
VTKM_CONT
MeshBoundary3D(vtkm::Id nrows,
vtkm::Id ncols,
vtkm::Id nslices,
const IdArrayType& sortOrder,
vtkm::cont::Token& token)
: MeshStructure(mesh_dem::MeshStructure3D<DeviceTag>(nrows, ncols, nslices))
{
this->SortOrderPortal = sortOrder.PrepareForInput(DeviceTag(), token);
}
VTKM_EXEC_CONT
bool liesOnBoundary(const vtkm::Id index) const
{
vtkm::Id meshSortOrderValue = this->SortOrderPortal.Get(index);
const vtkm::Id row = this->MeshStructure.VertexRow(meshSortOrderValue);
const vtkm::Id col = this->MeshStructure.VertexColumn(meshSortOrderValue);
const vtkm::Id sli = this->MeshStructure.VertexSlice(meshSortOrderValue);
return (row == 0) || (col == 0) || (sli == 0) || (row == this->MeshStructure.NumRows - 1) ||
(col == this->MeshStructure.NumColumns - 1) || (sli == this->MeshStructure.NumSlices - 1);
}
protected:
// 3D Mesh size parameters
mesh_dem::MeshStructure3D<DeviceTag> MeshStructure;
SortOrderPortalType SortOrderPortal;
};
class MeshBoundary3DExec : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_EXEC_CONT
MeshBoundary3DExec(vtkm::Id nrows,
vtkm::Id ncols,
vtkm::Id nslices,
const IdArrayType& inSortOrder)
: NumRows(nrows)
, NumColumns(ncols)
, NumSlices(nslices)
, SortOrder(inSortOrder)
{
}
VTKM_CONT
template <typename DeviceTag>
MeshBoundary3D<DeviceTag> PrepareForExecution(DeviceTag, vtkm::cont::Token& token) const
{
return MeshBoundary3D<DeviceTag>(
this->NumRows, this->NumColumns, this->NumSlices, this->SortOrder, token);
}
protected:
// 3D Mesh size parameters
vtkm::Id NumRows;
vtkm::Id NumColumns;
vtkm::Id NumSlices;
const IdArrayType& SortOrder;
};
} // namespace contourtree_augmented
} // worklet
} // vtkm
#endif

6
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/CMakeLists.txt → vtkm/worklet/contourtree_augmented/meshtypes/CMakeLists.txt
View File

@ -8,9 +8,9 @@
## PURPOSE. See the above copyright notice for more information.
##============================================================================
set(headers
Freudenthal_2D_Triangulation.h
Freudenthal_3D_Triangulation.h
MarchingCubes_3D_Triangulation.h
DataSetMeshTriangulation2DFreudenthal.h
DataSetMeshTriangulation3DFreudenthal.h
DataSetMeshTriangulation3DMarchingCubes.h
ContourTreeMesh.h
MeshStructureFreudenthal2D.h
MeshStructureFreudenthal3D.h

170
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/ContourTreeMesh.h → vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h
View File

@ -74,22 +74,23 @@
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/contourtree_augmented/ArrayTransforms.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ArcComparator.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedOtherStartIndexNNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedSimulatedSimplicityIndexComparator.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVector.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVectorDifferentFromNext.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CompressNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ComputeMaxNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/FindStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/MergeCombinedOtherStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/SubtractAssignWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/UpdateCombinedNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h> // This is needed only as an unused default argument.
#include <vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ArcComparator.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedOtherStartIndexNNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedSimulatedSimplicityIndexComparator.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVector.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVectorDifferentFromNext.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CompressNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ComputeMaxNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeCombinedOtherStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/SubtractAssignWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/UpdateCombinedNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h> // TODO remove should not be needed
@ -143,7 +144,7 @@ public:
const vtkm::cont::ArrayHandle<FieldType>& values,
const IdArrayType& inGlobalMeshIndex);
// Construct a ContourTreeMesh from nodes/arcs and another ContourTreeMesh (instead of a Mesh_DEM_Triangulation)
// Construct a ContourTreeMesh from nodes/arcs and another ContourTreeMesh (instead of a DataSetMesh)
// nodes/arcs: From the contour tree
// ContourTreeMesh: the contour tree mesh used to compute the contour tree described by nodes/arcs
ContourTreeMesh(const IdArrayType& nodes,
@ -185,8 +186,8 @@ public:
static const int MAX_OUTDEGREE = 20;
vtkm::Id NumVertices;
// TODO we should be able to remove this one, but we need to figure out what we need to return in the worklet instead
IdArrayType SortOrder;
vtkm::cont::ArrayHandleIndex SortOrder;
vtkm::cont::ArrayHandleIndex SortIndices;
vtkm::cont::ArrayHandle<FieldType> SortedValues;
IdArrayType GlobalMeshIndex;
// neighbours stores for each vertex the indices of its neighbours. For each vertex
@ -202,12 +203,14 @@ public:
// the maximum number of neighbours of a vertex
vtkm::Id MaxNeighbours;
// Print Contents
void PrintContent(std::ostream& outStream = std::cout) const;
// Debug print routine
void DebugPrint(const char* message, const char* fileName, long lineNum);
void DebugPrint(const char* message, const char* fileName, long lineNum) const;
// Get boundary execution object
MeshBoundaryContourTreeMeshExec GetMeshBoundaryExecutionObject(vtkm::Id totalNRows,
vtkm::Id totalNCols,
MeshBoundaryContourTreeMeshExec GetMeshBoundaryExecutionObject(vtkm::Id3 globalSize,
vtkm::Id3 minIdx,
vtkm::Id3 maxIdx) const;
@ -216,6 +219,42 @@ public:
MeshBoundaryContourTreeMeshExec* meshBoundaryExecObj //input
) const;
/// copies the global IDs for a set of sort IDs
/// notice that the sort ID is the same as the mesh ID for the ContourTreeMesh class.
/// To reduce memory usage we here use a fancy array handle rather than copy data
/// as is needed for the DataSetMesh types.
/// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
/// @param[in] sortIds Array with sort Ids to be converted from local to global Ids
/// @param[in] localToGlobalIdRelabeler This parameter is here only for
/// consistency with the DataSetMesh types but is not
/// used here and as such can simply be set to nullptr
inline vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> GetGlobalIdsFromSortIndices(
const IdArrayType& sortIds,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler =
nullptr) const
{ // GetGlobalIDsFromSortIndices()
(void)localToGlobalIdRelabeler; // avoid compiler warning
return vtkm::cont::make_ArrayHandlePermutation(sortIds, this->GlobalMeshIndex);
} // GetGlobalIDsFromSortIndices()
/// copies the global IDs for a set of mesh IDs
/// notice that the sort ID is the same as the mesh ID for the ContourTreeMesh class.
/// To reduce memory usage we here use a fancy array handle rather than copy data
/// as is needed for the DataSetMesh types.
/// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
/// @param[in] meshIds Array with mesh Ids to be converted from local to global Ids
/// @param[in] localToGlobalIdRelabeler This parameter is here only for
/// consistency with the DataSetMesh types but is not
/// used here and as such can simply be set to nullptr
inline vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> GetGlobalIdsFromMeshIndices(
const IdArrayType& meshIds,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler =
nullptr) const
{ // GetGlobalIDsFromMeshIndices()
(void)localToGlobalIdRelabeler; // avoid compiler warning
return vtkm::cont::make_ArrayHandlePermutation(meshIds, this->GlobalMeshIndex);
} // GetGlobalIDsFromMeshIndices()
private:
vtkm::cont::Invoker Invoke;
@ -237,10 +276,25 @@ private:
}; // ContourTreeMesh
// print content
template <typename FieldType>
void ContourTreeMesh<FieldType>::PrintContent(std::ostream& outStream /*= std::cout*/) const
{ // PrintContent()
PrintHeader(this->NumVertices, outStream);
//PrintIndices("SortOrder", SortOrder, outStream);
PrintValues("SortedValues", SortedValues, -1, outStream);
PrintIndices("GlobalMeshIndex", GlobalMeshIndex, -1, outStream);
PrintIndices("Neighbours", Neighbours, -1, outStream);
PrintIndices("FirstNeighbour", FirstNeighbour, -1, outStream);
outStream << "MaxNeighbours=" << MaxNeighbours << std::endl;
outStream << "mGetMax=" << mGetMax << std::endl;
} // DebugPrint()
// debug routine
template <typename FieldType>
void ContourTreeMesh<FieldType>::DebugPrint(const char* message, const char* fileName, long lineNum)
void ContourTreeMesh<FieldType>::DebugPrint(const char* message,
const char* fileName,
long lineNum) const
{ // DebugPrint()
#ifdef DEBUG_PRINT
std::cout << "---------------------------" << std::endl;
@ -251,15 +305,7 @@ void ContourTreeMesh<FieldType>::DebugPrint(const char* message, const char* fil
std::cout << "---------------------------" << std::endl;
std::cout << std::endl;
PrintHeader(this->NumVertices);
PrintIndices("SortOrder", SortOrder);
PrintValues("SortedValues", SortedValues);
PrintIndices("GlobalMeshIndex", GlobalMeshIndex);
PrintIndices("Neighbours", Neighbours);
PrintIndices("FirstNeighbour", FirstNeighbour);
std::cout << "MaxNeighbours=" << MaxNeighbours << std::endl;
std::cout << "mGetMax=" << mGetMax << std::endl;
PrintContent(std::cout);
#else
(void)message;
(void)fileName;
@ -267,24 +313,24 @@ void ContourTreeMesh<FieldType>::DebugPrint(const char* message, const char* fil
#endif
} // DebugPrint()
// create the contour tree mesh from contour tree data
template <typename FieldType>
ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& arcs,
const IdArrayType& inSortOrder,
const vtkm::cont::ArrayHandle<FieldType>& values,
const IdArrayType& inGlobalMeshIndex)
: SortOrder(inSortOrder)
: SortOrder()
, SortedValues()
, GlobalMeshIndex(inGlobalMeshIndex)
, Neighbours()
, FirstNeighbour()
{
this->NumVertices = this->SortOrder.GetNumberOfValues();
this->NumVertices = inSortOrder.GetNumberOfValues();
// Initalize the SortedIndices as a smart array handle
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
// values permuted by SortOrder to sort the values
auto permutedValues = vtkm::cont::make_ArrayHandlePermutation(this->SortOrder, values);
auto permutedValues = vtkm::cont::make_ArrayHandlePermutation(inSortOrder, values);
// TODO check if we actually need to make this copy here. we could just store the permutedValues array to save memory
vtkm::cont::Algorithm::Copy(permutedValues, this->SortedValues);
this->InitialiseNeighboursFromArcs(arcs);
@ -311,11 +357,10 @@ ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
inSortOrder);
auto permutedValues = vtkm::cont::make_ArrayHandlePermutation(permutedSortOrder, values);
vtkm::cont::Algorithm::Copy(permutedValues, this->SortedValues);
vtkm::cont::Algorithm::Copy(
permutedSortOrder,
this
->SortOrder); // TODO Check if the SortOrder needs to be set form the input or the permutted SortOrder
// Initalize the SortedIndices as a smart array handle
this->NumVertices = this->SortedValues.GetNumberOfValues();
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->InitialiseNeighboursFromArcs(arcs);
#ifdef DEBUG_PRINT
// Print the contents fo this for debugging
@ -326,13 +371,15 @@ ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
template <typename FieldType>
ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& arcs,
const ContourTreeMesh<FieldType>& mesh)
: SortOrder(mesh.SortOrder)
, SortedValues(mesh.SortedValues)
: SortedValues(mesh.SortedValues)
, GlobalMeshIndex(mesh.GlobalMeshIndex)
, Neighbours()
, FirstNeighbour()
{
// Initalize the SortedIndices as a smart array handle
this->NumVertices = this->SortedValues.GetNumberOfValues();
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->InitialiseNeighboursFromArcs(arcs);
#ifdef DEBUG_PRINT
// Print the contents fo this for debugging
@ -345,8 +392,7 @@ template <typename FieldType>
ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
const IdArrayType& arcs,
const ContourTreeMesh<FieldType>& mesh)
: SortOrder(mesh.SortOrder)
, Neighbours()
: Neighbours()
, FirstNeighbour()
{
// Initatlize the global mesh index with the GlobalMeshIndex permutted by the nodes
@ -358,6 +404,8 @@ ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
vtkm::cont::Algorithm::Copy(permutedSortedValues, this->SortedValues);
// Initialize the neighbours from the arcs
this->NumVertices = this->SortedValues.GetNumberOfValues();
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->InitialiseNeighboursFromArcs(arcs);
#ifdef DEBUG_PRINT
// Print the contents fo this for debugging
@ -479,7 +527,7 @@ struct NotNoSuchElement
template <typename FieldType>
template <typename DeviceTag>
void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& other)
{
{ // Merge With
#ifdef DEBUG_PRINT
this->DebugPrint("THIS ContourTreeMesh", __FILE__, __LINE__);
other.DebugPrint("OTHER ContourTreeMesh", __FILE__, __LINE__);
@ -495,8 +543,8 @@ void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& other)
//auto allGlobalIndices = CombinedVector<FieldType(this->thisGlobalMeshIndex, other.GlobalMeshIndex);
// Create combined sort order
IdArrayType
overallSortOrder; // TODO This vector could potentially be implemented purely as a smart array handle to reduce memory usage
// TODO This vector could potentially be implemented purely as a smart array handle to reduce memory usage
IdArrayType overallSortOrder;
overallSortOrder.Allocate(this->NumVertices + other.NumVertices);
{ // Create a new scope so that the following two vectors get deleted when leaving the scope
@ -725,7 +773,8 @@ void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& other)
this->Neighbours = combinedNeighbours;
this->FirstNeighbour = combinedFirstNeighbour;
this->NumVertices = SortedValues.GetNumberOfValues();
// TODO Do we need to set the sort order as well?
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
// Re-compute maximum number of neigbours
ComputeMaxNeighbours();
@ -734,14 +783,13 @@ void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& other)
// Print the contents fo this for debugging
DebugPrint("ContourTreeMeshes merged", __FILE__, __LINE__);
#endif
}
} // Merge With
template <typename FieldType>
void ContourTreeMesh<FieldType>::Save(const char* filename) const
{
std::ofstream os(filename);
SaveVector(os, this->SortOrder);
SaveVector(os, this->SortedValues);
SaveVector(os, this->GlobalMeshIndex);
SaveVector(os, this->Neighbours);
@ -752,12 +800,14 @@ template <typename FieldType>
void ContourTreeMesh<FieldType>::Load(const char* filename)
{
std::ifstream is(filename);
LoadVector(is, this->SortOrder);
LoadVector(is, this->SortedValues);
LoadVector(is, this->GlobalMeshIndex);
LoadVector(is, this->Neighbours);
LoadVector(is, this->FirstNeighbour);
this->ComputeMaxNeighbours();
this->NumVertices = this->SortedValues.GetNumberOfValues();
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
}
template <typename FieldType>
@ -766,10 +816,13 @@ void ContourTreeMesh<FieldType>::SaveVector(std::ostream& os,
const vtkm::cont::ArrayHandle<ValueType>& vec) const
{
vtkm::Id numVals = vec.GetNumberOfValues();
os.write(reinterpret_cast<const char*>(&numVals), sizeof(ValueType));
//os.write(reinterpret_cast<const char*>(&numVals), sizeof(ValueType));
os << numVals << ": ";
auto vecPortal = vec.ReadPortal();
for (vtkm::Id i = 0; i < numVals; ++i)
os.write(reinterpret_cast<const char*>(vecPortal.Get(i)), sizeof(ValueType));
os << vecPortal.Get(i) << " ";
//os.write(reinterpret_cast<const char*>(vecPortal.Get(i)), sizeof(ValueType));
os << std::endl;
}
template <typename FieldType>
@ -791,13 +844,11 @@ void ContourTreeMesh<FieldType>::LoadVector(std::istream& is,
template <typename FieldType>
MeshBoundaryContourTreeMeshExec ContourTreeMesh<FieldType>::GetMeshBoundaryExecutionObject(
vtkm::Id totalNRows,
vtkm::Id totalNCols,
vtkm::Id3 globalSize,
vtkm::Id3 minIdx,
vtkm::Id3 maxIdx) const
{
return MeshBoundaryContourTreeMeshExec(
this->GlobalMeshIndex, totalNRows, totalNCols, minIdx, maxIdx);
return MeshBoundaryContourTreeMeshExec(this->GlobalMeshIndex, globalSize, minIdx, maxIdx);
}
template <typename FieldType>
@ -824,7 +875,6 @@ void ContourTreeMesh<FieldType>::GetBoundaryVertices(
vtkm::cont::Algorithm::Copy(boundaryVertexArray, boundarySortIndexArray);
}
} // namespace contourtree_augmented
} // worklet
} // vtkm

64
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/Freudenthal_2D_Triangulation.h → vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation2DFreudenthal.h
View File

@ -50,16 +50,16 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_2d_freudenthal_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_2d_freudenthal_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_2d_freudenthal_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_2d_freudenthal_h
#include <cstdlib>
#include <vtkm/Types.h>
#include <vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureFreudenthal2D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundary2D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary2D.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureFreudenthal2D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundary2D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary2D.h>
#include <vtkm/cont/ExecutionObjectBase.h>
@ -70,14 +70,14 @@ namespace worklet
namespace contourtree_augmented
{
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation_2D_Freudenthal
: public Mesh_DEM_Triangulation_2D<T, StorageType>
class DataSetMeshTriangulation2DFreudenthal
: public DataSetMesh
, public vtkm::cont::ExecutionObjectBase
{ // class Mesh_DEM_Triangulation
{ // class DataSetMeshTriangulation
public:
// Constants and case tables
m2d_freudenthal::EdgeBoundaryDetectionMasksType EdgeBoundaryDetectionMasks;
static constexpr int MAX_OUTDEGREE = 3;
//Mesh dependent helper functions
void SetPrepareForExecutionBehavior(bool getMax);
@ -86,7 +86,7 @@ public:
MeshStructureFreudenthal2D<DeviceTag> PrepareForExecution(DeviceTag,
vtkm::cont::Token& token) const;
Mesh_DEM_Triangulation_2D_Freudenthal(vtkm::Id ncols, vtkm::Id nrows);
DataSetMeshTriangulation2DFreudenthal(vtkm::Id2 meshSize);
MeshBoundary2DExec GetMeshBoundaryExecutionObject() const;
@ -98,39 +98,30 @@ public:
private:
bool UseGetMax; // Define the behavior ofr the PrepareForExecution function
}; // class Mesh_DEM_Triangulation
}; // class DataSetMeshTriangulation
// creates input mesh
template <typename T, typename StorageType>
Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::Mesh_DEM_Triangulation_2D_Freudenthal(
vtkm::Id ncols,
vtkm::Id nrows)
: Mesh_DEM_Triangulation_2D<T, StorageType>(ncols, nrows)
DataSetMeshTriangulation2DFreudenthal::DataSetMeshTriangulation2DFreudenthal(vtkm::Id2 meshSize)
: DataSetMesh(vtkm::Id3{ meshSize[0], meshSize[1], 1 })
, EdgeBoundaryDetectionMasks{ vtkm::cont::make_ArrayHandle(
m2d_freudenthal::EdgeBoundaryDetectionMasks,
m2d_freudenthal::N_INCIDENT_EDGES,
vtkm::CopyFlag::Off) }
{
this->EdgeBoundaryDetectionMasks =
vtkm::cont::make_ArrayHandle(m2d_freudenthal::EdgeBoundaryDetectionMasks,
m2d_freudenthal::N_INCIDENT_EDGES,
vtkm::CopyFlag::Off);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::SetPrepareForExecutionBehavior(
bool getMax)
void DataSetMeshTriangulation2DFreudenthal::SetPrepareForExecutionBehavior(bool getMax)
{
this->UseGetMax = getMax;
}
// Get VTKM execution object that represents the structure of the mesh and provides the mesh helper functions on the device
template <typename T, typename StorageType>
template <typename DeviceTag>
MeshStructureFreudenthal2D<DeviceTag>
Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::PrepareForExecution(
MeshStructureFreudenthal2D<DeviceTag> DataSetMeshTriangulation2DFreudenthal::PrepareForExecution(
DeviceTag,
vtkm::cont::Token& token) const
{
return MeshStructureFreudenthal2D<DeviceTag>(this->NumColumns,
this->NumRows,
return MeshStructureFreudenthal2D<DeviceTag>(vtkm::Id2{ this->MeshSize[0], this->MeshSize[1] },
m2d_freudenthal::N_INCIDENT_EDGES,
this->UseGetMax,
this->SortIndices,
@ -139,23 +130,19 @@ Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::PrepareForExecution(
token);
}
template <typename T, typename StorageType>
MeshBoundary2DExec
Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::GetMeshBoundaryExecutionObject() const
MeshBoundary2DExec DataSetMeshTriangulation2DFreudenthal::GetMeshBoundaryExecutionObject() const
{
return MeshBoundary2DExec(this->NumColumns, this->NumRows, this->SortOrder);
return MeshBoundary2DExec(vtkm::Id2{ this->MeshSize[0], this->MeshSize[1] }, this->SortIndices);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::GetBoundaryVertices(
void DataSetMeshTriangulation2DFreudenthal::GetBoundaryVertices(
IdArrayType& boundaryVertexArray, // output
IdArrayType& boundarySortIndexArray, // output
MeshBoundary2DExec*
meshBoundaryExecObj // optional input, included for consistency with ContourTreeMesh
) const
{
vtkm::Id numBoundary = 2 * this->NumRows + 2 * this->NumColumns - 4;
vtkm::Id numBoundary = 2 * this->MeshSize[1] + 2 * this->MeshSize[0] - 4;
auto boundaryId = vtkm::cont::ArrayHandleIndex(numBoundary);
ComputeMeshBoundary2D computeMeshBoundary2dWorklet;
this->Invoke(computeMeshBoundary2dWorklet,
@ -168,7 +155,6 @@ void Mesh_DEM_Triangulation_2D_Freudenthal<T, StorageType>::GetBoundaryVertices(
);
}
} // namespace contourtree_augmented
} // worklet
} // vtkm

63
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/Freudenthal_3D_Triangulation.h → vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DFreudenthal.h
View File

@ -51,8 +51,8 @@
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_3d_freudenthal_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_3d_freudenthal_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_3d_freudenthal_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_3d_freudenthal_h
#include <cstdlib>
#include <vtkm/Types.h>
@ -60,11 +60,11 @@
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureFreudenthal3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/freudenthal_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureFreudenthal3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/freudenthal_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary3D.h>
namespace vtkm
{
@ -73,16 +73,16 @@ namespace worklet
namespace contourtree_augmented
{
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation_3D_Freudenthal
: public Mesh_DEM_Triangulation_3D<T, StorageType>
class DataSetMeshTriangulation3DFreudenthal
: public DataSetMesh
, public vtkm::cont::ExecutionObjectBase
{ // class Mesh_DEM_Triangulation
{ // class DataSetMeshTriangulation3DFreudenthal
public:
// Constants and case tables
m3d_freudenthal::EdgeBoundaryDetectionMasksType EdgeBoundaryDetectionMasks;
m3d_freudenthal::NeighbourOffsetsType NeighbourOffsets;
m3d_freudenthal::LinkComponentCaseTableType LinkComponentCaseTable;
static constexpr int MAX_OUTDEGREE = 6; // True for Freudenthal and Marching Cubes
// Mesh helper functions
void SetPrepareForExecutionBehavior(bool getMax);
@ -91,7 +91,7 @@ public:
MeshStructureFreudenthal3D<DeviceTag> PrepareForExecution(DeviceTag,
vtkm::cont::Token& token) const;
Mesh_DEM_Triangulation_3D_Freudenthal(vtkm::Id ncols, vtkm::Id nrows, vtkm::Id nslices);
DataSetMeshTriangulation3DFreudenthal(vtkm::Id3 meshSize);
MeshBoundary3DExec GetMeshBoundaryExecutionObject() const;
@ -103,15 +103,11 @@ public:
private:
bool UseGetMax; // Define the behavior ofr the PrepareForExecution function
}; // class Mesh_DEM_Triangulation
}; // class DataSetMeshTriangulation
// creates input mesh
template <typename T, typename StorageType>
Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::Mesh_DEM_Triangulation_3D_Freudenthal(
vtkm::Id ncols,
vtkm::Id nrows,
vtkm::Id nslices)
: Mesh_DEM_Triangulation_3D<T, StorageType>(ncols, nrows, nslices)
DataSetMeshTriangulation3DFreudenthal::DataSetMeshTriangulation3DFreudenthal(vtkm::Id3 meshSize)
: DataSetMesh(meshSize)
{
// Initialize the case tables in vtkm
@ -127,25 +123,18 @@ Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::Mesh_DEM_Triangulation_3D
vtkm::CopyFlag::Off);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::SetPrepareForExecutionBehavior(
bool getMax)
void DataSetMeshTriangulation3DFreudenthal::SetPrepareForExecutionBehavior(bool getMax)
{
this->UseGetMax = getMax;
}
// Get VTKM execution object that represents the structure of the mesh and provides the mesh helper functions on the device
template <typename T, typename StorageType>
template <typename DeviceTag>
MeshStructureFreudenthal3D<DeviceTag>
Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::PrepareForExecution(
MeshStructureFreudenthal3D<DeviceTag> DataSetMeshTriangulation3DFreudenthal::PrepareForExecution(
DeviceTag,
vtkm::cont::Token& token) const
{
return MeshStructureFreudenthal3D<DeviceTag>(this->NumColumns,
this->NumRows,
this->NumSlices,
return MeshStructureFreudenthal3D<DeviceTag>(this->MeshSize,
m3d_freudenthal::N_INCIDENT_EDGES,
this->UseGetMax,
this->SortIndices,
@ -156,24 +145,20 @@ Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::PrepareForExecution(
token);
}
template <typename T, typename StorageType>
MeshBoundary3DExec
Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::GetMeshBoundaryExecutionObject() const
MeshBoundary3DExec DataSetMeshTriangulation3DFreudenthal::GetMeshBoundaryExecutionObject() const
{
return MeshBoundary3DExec(this->NumColumns, this->NumRows, this->NumSlices, this->SortOrder);
return MeshBoundary3DExec(this->MeshSize, this->SortIndices);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D_Freudenthal<T, StorageType>::GetBoundaryVertices(
void DataSetMeshTriangulation3DFreudenthal::GetBoundaryVertices(
IdArrayType& boundaryVertexArray, // output
IdArrayType& boundarySortIndexArray, // output
MeshBoundary3DExec* meshBoundaryExecObj // input
) const
{
vtkm::Id numBoundary = 2 * this->NumRows * this->NumColumns // xy faces
+ 2 * this->NumRows * (this->NumSlices - 2) // yz faces - excluding vertices on xy
+ 2 * (this->NumColumns - 2) * (this->NumSlices - 2); // xz face interiors
vtkm::Id numBoundary = 2 * this->MeshSize[1] * this->MeshSize[0] // xy faces
+ 2 * this->MeshSize[1] * (this->MeshSize[2] - 2) // yz faces - excluding vertices on xy
+ 2 * (this->MeshSize[0] - 2) * (this->MeshSize[2] - 2); // xz face interiors
auto boundaryId = vtkm::cont::ArrayHandleIndex(numBoundary);
ComputeMeshBoundary3D computeMeshBoundary3dWorklet;
this->Invoke(computeMeshBoundary3dWorklet,

70
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MarchingCubes_3D_Triangulation.h → vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DMarchingCubes.h
View File

@ -51,8 +51,8 @@
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_3d_marchingcubes_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_triangulation_3d_marchingcubes_h
#ifndef vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_3d_marchingcubes_h
#define vtk_m_worklet_contourtree_augmented_data_set_mesh_triangulation_3d_marchingcubes_h
#include <cstdlib>
#include <vtkm/Types.h>
@ -60,11 +60,11 @@
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/contourtree_augmented/Mesh_DEM_Triangulation.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureMarchingCubes.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/marchingcubes_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureMarchingCubes.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/marchingcubes_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundary3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary3D.h>
namespace vtkm
{
@ -73,11 +73,10 @@ namespace worklet
namespace contourtree_augmented
{
template <typename T, typename StorageType>
class Mesh_DEM_Triangulation_3D_MarchingCubes
: public Mesh_DEM_Triangulation_3D<T, StorageType>
class DataSetMeshTriangulation3DMarchingCubes
: public DataSetMesh
, public vtkm::cont::ExecutionObjectBase
{ // class Mesh_DEM_Triangulation
{ // class DataSetMeshTriangulation3DMarchingCubes
public:
//Constants and case tables
@ -87,6 +86,7 @@ public:
m3d_marchingcubes::LinkVertexConnectionsType LinkVertexConnectionsEighteen;
m3d_marchingcubes::InCubeConnectionsType InCubeConnectionsSix;
m3d_marchingcubes::InCubeConnectionsType InCubeConnectionsEighteen;
static constexpr int MAX_OUTDEGREE = 6; // True for Freudenthal and Marching Cubes
// mesh depended helper functions
void SetPrepareForExecutionBehavior(bool getMax);
@ -95,27 +95,24 @@ public:
MeshStructureMarchingCubes<DeviceTag> PrepareForExecution(DeviceTag,
vtkm::cont::Token& token) const;
Mesh_DEM_Triangulation_3D_MarchingCubes(vtkm::Id ncols, vtkm::Id nrows, vtkm::Id nslices);
DataSetMeshTriangulation3DMarchingCubes(vtkm::Id3 meshSize);
MeshBoundary3DExec GetMeshBoundaryExecutionObject() const;
void GetBoundaryVertices(IdArrayType& boundaryVertexArray, // output
IdArrayType& boundarySortIndexArray, // output
MeshBoundary3DExec* meshBoundaryExecObj =
NULL // optional input, included for consistency with ContourTreeMesh
void GetBoundaryVertices(
IdArrayType& boundaryVertexArray, // output
IdArrayType& boundarySortIndexArray, // output
MeshBoundary3DExec* meshBoundaryExecObj =
nullptr // optional input, included for consistency with ContourTreeMesh
) const;
private:
bool UseGetMax; // Define the behavior ofr the PrepareForExecution function
}; // class Mesh_DEM_Triangulation
}; // class DataSetMesh_Triangulation
// creates input mesh
template <typename T, typename StorageType>
Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::Mesh_DEM_Triangulation_3D_MarchingCubes(
vtkm::Id ncols,
vtkm::Id nrows,
vtkm::Id nslices)
: Mesh_DEM_Triangulation_3D<T, StorageType>(ncols, nrows, nslices)
DataSetMeshTriangulation3DMarchingCubes::DataSetMeshTriangulation3DMarchingCubes(vtkm::Id3 meshSize)
: DataSetMesh(meshSize)
{
// Initialize the case tables in vtkm
@ -157,24 +154,18 @@ Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::Mesh_DEM_Triangulation_
vtkm::CopyFlag::Off);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::SetPrepareForExecutionBehavior(
bool getMax)
void DataSetMeshTriangulation3DMarchingCubes::SetPrepareForExecutionBehavior(bool getMax)
{
this->UseGetMax = getMax;
}
// Get VTKM execution object that represents the structure of the mesh and provides the mesh helper functions on the device
template <typename T, typename StorageType>
template <typename DeviceTag>
MeshStructureMarchingCubes<DeviceTag>
Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::PrepareForExecution(
MeshStructureMarchingCubes<DeviceTag> DataSetMeshTriangulation3DMarchingCubes::PrepareForExecution(
DeviceTag,
vtkm::cont::Token& token) const
{
return MeshStructureMarchingCubes<DeviceTag>(this->NumColumns,
this->NumRows,
this->NumSlices,
return MeshStructureMarchingCubes<DeviceTag>(this->MeshSize,
this->UseGetMax,
this->SortIndices,
this->SortOrder,
@ -187,23 +178,20 @@ Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::PrepareForExecution(
token);
}
template <typename T, typename StorageType>
MeshBoundary3DExec
Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::GetMeshBoundaryExecutionObject() const
MeshBoundary3DExec DataSetMeshTriangulation3DMarchingCubes::GetMeshBoundaryExecutionObject() const
{
return MeshBoundary3DExec(this->NumColumns, this->NumRows, this->NumSlices, this->SortOrder);
return MeshBoundary3DExec(this->MeshSize, this->SortOrder);
}
template <typename T, typename StorageType>
void Mesh_DEM_Triangulation_3D_MarchingCubes<T, StorageType>::GetBoundaryVertices(
void DataSetMeshTriangulation3DMarchingCubes::GetBoundaryVertices(
IdArrayType& boundaryVertexArray, // output
IdArrayType& boundarySortIndexArray, // output
MeshBoundary3DExec* meshBoundaryExecObj // input
) const
{
vtkm::Id numBoundary = 2 * this->NumRows * this->NumColumns // xy faces
+ 2 * this->NumRows * (this->NumSlices - 2) // yz faces - excluding vertices on xy
+ 2 * (this->NumColumns - 2) * (this->NumSlices - 2); // xz face interiors
vtkm::Id numBoundary = 2 * this->MeshSize[1] * this->MeshSize[0] // xy faces
+ 2 * this->MeshSize[1] * (this->MeshSize[2] - 2) // yz faces - excluding vertices on xy
+ 2 * (this->MeshSize[0] - 2) * (this->MeshSize[2] - 2); // xz face interiors
auto boundaryId = vtkm::cont::ArrayHandleIndex(numBoundary);
ComputeMeshBoundary3D computeMeshBoundary3dWorklet;
this->Invoke(computeMeshBoundary3dWorklet,

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureContourTreeMesh.h → vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureContourTreeMesh.h
View File

47
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureFreudenthal2D.h → vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureFreudenthal2D.h
View File

@ -50,15 +50,15 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureFreudenthal2D_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureFreudenthal2D_h
#ifndef vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureFreudenthal2D_h
#define vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureFreudenthal2D_h
#include <vtkm/Pair.h>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure2D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/freudenthal_2D/Types.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure2D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/freudenthal_2D/Types.h>
namespace vtkm
@ -70,7 +70,7 @@ namespace contourtree_augmented
// Worklet for computing the sort indices from the sort order
template <typename DeviceAdapter>
class MeshStructureFreudenthal2D : public mesh_dem::MeshStructure2D<DeviceAdapter>
class MeshStructureFreudenthal2D : public data_set_mesh::MeshStructure2D<DeviceAdapter>
{
public:
using SortIndicesPortalType =
@ -82,7 +82,7 @@ public:
// Default constucture. Needed for the CUDA built to work
VTKM_EXEC_CONT
MeshStructureFreudenthal2D()
: mesh_dem::MeshStructure2D<DeviceAdapter>()
: data_set_mesh::MeshStructure2D<DeviceAdapter>()
, GetMax(false)
, NumIncidentEdges(m2d_freudenthal::N_INCIDENT_EDGES)
{
@ -90,15 +90,14 @@ public:
// Main constructor used in the code
MeshStructureFreudenthal2D(
vtkm::Id ncols,
vtkm::Id nrows,
vtkm::Id2 meshSize,
vtkm::Int32 nincident_edges,
bool getmax,
const IdArrayType& sortIndices,
const IdArrayType& SortOrder,
const m2d_freudenthal::EdgeBoundaryDetectionMasksType& EdgeBoundaryDetectionMasksIn,
vtkm::cont::Token& token)
: mesh_dem::MeshStructure2D<DeviceAdapter>(ncols, nrows)
: data_set_mesh::MeshStructure2D<DeviceAdapter>(meshSize)
, GetMax(getmax)
, NumIncidentEdges(nincident_edges)
{
@ -118,17 +117,17 @@ public:
switch (edgeNo)
{
case 0:
return this->SortIndicesPortal.Get(meshIndex + 1); // row , col + 1
return this->SortIndicesPortal.Get(meshIndex + 1); // [1] , [0] + 1
case 1:
return this->SortIndicesPortal.Get(meshIndex + this->NumColumns + 1); // row + 1, col + 1
return this->SortIndicesPortal.Get(meshIndex + this->MeshSize[0] + 1); // [1] + 1, [0] + 1
case 2:
return this->SortIndicesPortal.Get(meshIndex + this->NumColumns); // row + 1, col
return this->SortIndicesPortal.Get(meshIndex + this->MeshSize[0]); // [1] + 1, [0]
case 3:
return this->SortIndicesPortal.Get(meshIndex - 1); // row , col - 1
return this->SortIndicesPortal.Get(meshIndex - 1); // [1] , [0] - 1
case 4:
return this->SortIndicesPortal.Get(meshIndex - this->NumColumns - 1); // row - 1, col - 1
return this->SortIndicesPortal.Get(meshIndex - this->MeshSize[0] - 1); // [1] - 1, [0] - 1
case 5:
return this->SortIndicesPortal.Get(meshIndex - this->NumColumns); // row - 1, col
return this->SortIndicesPortal.Get(meshIndex - this->MeshSize[0]); // [1] - 1, [0]
default:
return -1; // TODO How to generate a meaningful error message from a device (in particular when using CUDA?)
}
@ -154,11 +153,10 @@ public:
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
// get the row and column
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id2 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0);
// in what follows, the boundary conditions always reset wasAscent
for (vtkm::Id edgeNo = 0; edgeNo < this->NumIncidentEdges; edgeNo++)
@ -194,11 +192,10 @@ public:
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
// get the row and column
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id2 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0);
// and initialise the mask
vtkm::Id neighbourhoodMask = 0;

49
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureFreudenthal3D.h → vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureFreudenthal3D.h
View File

@ -50,16 +50,16 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureFreudenthal3D_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureFreudenthal3D_h
#ifndef vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureFreudenthal3D_h
#define vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureFreudenthal3D_h
#include <vtkm/Pair.h>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/freudenthal_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/freudenthal_3D/Types.h>
namespace vtkm
@ -71,7 +71,7 @@ namespace contourtree_augmented
// Worklet for computing the sort indices from the sort order
template <typename DeviceAdapter>
class MeshStructureFreudenthal3D : public mesh_dem::MeshStructure3D<DeviceAdapter>
class MeshStructureFreudenthal3D : public data_set_mesh::MeshStructure3D<DeviceAdapter>
{
public:
using SortIndicesPortalType =
@ -92,7 +92,7 @@ public:
// Default constructor needed to make the CUDA build work
VTKM_EXEC_CONT
MeshStructureFreudenthal3D()
: mesh_dem::MeshStructure3D<DeviceAdapter>()
: data_set_mesh::MeshStructure3D<DeviceAdapter>()
, GetMax(false)
, NumIncidentEdge(m3d_freudenthal::N_INCIDENT_EDGES)
{
@ -100,9 +100,7 @@ public:
// Main constructore used in the code
MeshStructureFreudenthal3D(
vtkm::Id ncols,
vtkm::Id nrows,
vtkm::Id nslices,
vtkm::Id3 meshSize,
vtkm::Id nincident_edges,
bool getmax,
const IdArrayType& sortIndices,
@ -111,7 +109,7 @@ public:
const m3d_freudenthal::NeighbourOffsetsType& neighbourOffsetsIn,
const m3d_freudenthal::LinkComponentCaseTableType& linkComponentCaseTableIn,
vtkm::cont::Token& token)
: mesh_dem::MeshStructure3D<DeviceAdapter>(ncols, nrows, nslices)
: data_set_mesh::MeshStructure3D<DeviceAdapter>(meshSize)
, GetMax(getmax)
, NumIncidentEdge(nincident_edges)
{
@ -132,10 +130,12 @@ public:
inline vtkm::Id GetNeighbourIndex(vtkm::Id sortIndex, vtkm::Id edgeNo) const
{ // GetNeighbourIndex
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
// NOTE: Offsets are stored in "reversed" zyx [2][1][0] order (remaining artifact from
// using slices, rows, columns instead of xyz/[0][1][2])
return SortIndicesPortal.Get(meshIndex +
(NeighbourOffsetsPortal.Get(edgeNo)[0] * this->NumRows +
(NeighbourOffsetsPortal.Get(edgeNo)[0] * this->MeshSize[1] +
NeighbourOffsetsPortal.Get(edgeNo)[1]) *
this->NumColumns +
this->MeshSize[0] +
NeighbourOffsetsPortal.Get(edgeNo)[2]);
} // GetNeighbourIndex
@ -159,14 +159,11 @@ public:
using namespace m3d_freudenthal;
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
vtkm::Id slice = this->VertexSlice(meshIndex);
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((slice == 0) ? FrontBit : 0) |
((slice == this->NumSlices - 1) ? BackBit : 0) | ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id3 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0) | ((pos[2] == 0) ? FrontBit : 0) |
((pos[2] == this->MeshSize[2] - 1) ? BackBit : 0);
// in what follows, the boundary conditions always reset wasAscent
// loop downwards so that we pick the same edges as previous versions
@ -199,13 +196,11 @@ public:
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
// get the row and column
vtkm::Id slice = this->VertexSlice(meshIndex);
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((slice == 0) ? FrontBit : 0) |
((slice == this->NumSlices - 1) ? BackBit : 0) | ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id3 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0) | ((pos[2] == 0) ? FrontBit : 0) |
((pos[2] == this->MeshSize[2] - 1) ? BackBit : 0);
// Initialize "union find"
vtkm::Id caseNo = 0;

169
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/MeshStructureMarchingCubes.h → vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureMarchingCubes.h
View File

@ -50,16 +50,16 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureMarchingCubes_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_MeshStructureMarchingCubes_h
#ifndef vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureMarchingCubes_h
#define vtk_m_worklet_contourtree_augmented_meshtypes_MeshStructureMarchingCubes_h
#include <vtkm/Pair.h>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure3D.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/marchingcubes_3D/Types.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure3D.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/marchingcubes_3D/Types.h>
namespace vtkm
{
@ -70,7 +70,7 @@ namespace contourtree_augmented
// Worklet for computing the sort indices from the sort order
template <typename DeviceAdapter>
class MeshStructureMarchingCubes : public mesh_dem::MeshStructure3D<DeviceAdapter>
class MeshStructureMarchingCubes : public data_set_mesh::MeshStructure3D<DeviceAdapter>
{
public:
// EdgeBoundaryDetectionMasks types
@ -100,16 +100,14 @@ public:
// Default constructor needed to make the CUDA build work
VTKM_EXEC_CONT
MeshStructureMarchingCubes()
: mesh_dem::MeshStructure3D<DeviceAdapter>()
: data_set_mesh::MeshStructure3D<DeviceAdapter>()
, GetMax(false)
{
}
// Main constructore used in the code
MeshStructureMarchingCubes(
vtkm::Id ncols,
vtkm::Id nrows,
vtkm::Id nslices,
vtkm::Id3 meshSize,
bool getmax,
const IdArrayType& sortIndices,
const IdArrayType& sortOrder,
@ -120,7 +118,7 @@ public:
const m3d_marchingcubes::InCubeConnectionsType& InCubeConnectionsSixIn,
const m3d_marchingcubes::InCubeConnectionsType& InCubeConnectionsEighteenIn,
vtkm::cont::Token& token)
: mesh_dem::MeshStructure3D<DeviceAdapter>(ncols, nrows, nslices)
: data_set_mesh::MeshStructure3D<DeviceAdapter>(meshSize)
, GetMax(getmax)
{
this->SortIndicesPortal = sortIndices.PrepareForInput(DeviceAdapter(), token);
@ -150,84 +148,69 @@ public:
{
using namespace m3d_marchingcubes;
vtkm::Id meshIndex = this->SortOrderPortal.Get(sortIndex);
const vtkm::Id3 strides{ 1, this->MeshSize[0], this->MeshSize[0] * this->MeshSize[1] };
// GetNeighbourIndex
switch (nbrNo)
{
// Edge connected neighbours
case 0:
return SortIndicesPortal.Get(meshIndex -
(this->NumRows * this->NumColumns)); // { -1, 0, 0 }
case 1:
return SortIndicesPortal.Get(meshIndex - this->NumColumns); // { 0, -1, 0 }
case 2:
return SortIndicesPortal.Get(meshIndex - 1); // { 0, 0, -1 }
case 3:
return SortIndicesPortal.Get(meshIndex + 1); // { 0, 0, 1 }
case 4:
return SortIndicesPortal.Get(meshIndex + this->NumColumns); // { 0, 1, 0 }
case 0: // { 0, 0, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2]);
case 1: // { 0, -1, 0 }
return SortIndicesPortal.Get(meshIndex - strides[1]);
case 2: // { -1, 0, 0 }
return SortIndicesPortal.Get(meshIndex - strides[0]);
case 3: // { 1, 0, 0 }
return SortIndicesPortal.Get(meshIndex + strides[0]);
case 4: // { 0, 1, 0 }
return SortIndicesPortal.Get(meshIndex + strides[1]);
case 5: // { 0, 0, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2]);
// Face connected neighbours
case 5:
return SortIndicesPortal.Get(meshIndex +
(this->NumRows * this->NumColumns)); // { 1, 0, 0 }
case 6:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) -
this->NumColumns); // { -1, -1, 0 }
case 7:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) -
1); // { -1, 0, -1 }
case 8:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) +
1); // { -1, 0, 1 }
case 9:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) +
this->NumColumns); // { -1, 1, 0 }
case 10:
return SortIndicesPortal.Get(meshIndex - this->NumColumns - 1); // { 0, -1, -1 }
case 11:
return SortIndicesPortal.Get(meshIndex - this->NumColumns + 1); // { 0, -1, 1 }
case 12:
return SortIndicesPortal.Get(meshIndex + this->NumColumns - 1); // { 0, 1, -1 }
case 13:
return SortIndicesPortal.Get(meshIndex + this->NumColumns + 1); // { 0, 1, 1 }
case 14:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) -
this->NumColumns); // { 1, -1, 0 }
case 15:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) -
1); // { 1, 0, -1 }
case 16:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) +
1); // { 1, 0, 1 }
case 17:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) +
this->NumColumns); // { 1, 1, 0 }
case 6: // { 0, -1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] - strides[1]);
case 7: // { -1, 0, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] - strides[0]);
case 8: // { 1, 0, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] + strides[0]);
case 9: // { 0, 1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] + strides[1]);
case 10: // { -1, -1, 0 }
return SortIndicesPortal.Get(meshIndex - strides[1] - strides[0]);
case 11: // { 1, -1, 0 }
return SortIndicesPortal.Get(meshIndex - strides[1] + strides[0]);
case 12: // { -1, 1, 0 }
return SortIndicesPortal.Get(meshIndex + strides[1] - strides[0]);
case 13: // { 1, 1, 0 }
return SortIndicesPortal.Get(meshIndex + strides[1] + strides[0]);
case 14: // { 0, -1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] - strides[1]);
case 15: // { -1, 0, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] - 1);
case 16: // { 1, 0, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] + 1);
case 17: // { 0, 1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] + strides[1]);
// Diagonal connected neighbours
case 18:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) -
this->NumColumns - 1); // { -1, -1, -1 }
case 19:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) -
this->NumColumns + 1); // { -1, -1, 1 }
case 20:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) +
this->NumColumns - 1); // { -1, 1, -1 }
case 21:
return SortIndicesPortal.Get(meshIndex - (this->NumRows * this->NumColumns) +
this->NumColumns + 1); // { -1, 1, 1 }
case 22:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) -
this->NumColumns - 1); // { 1, -1, -1 }
case 23:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) -
this->NumColumns + 1); // { 1, -1, 1 }
case 24:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) +
this->NumColumns - 1); // { 1, 1, -1 }
case 25:
return SortIndicesPortal.Get(meshIndex + (this->NumRows * this->NumColumns) +
this->NumColumns + 1); // { 1, 1, 1 }
case 18: // { -1, -1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] - strides[1] - strides[0]);
case 19: // { 1, -1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] - strides[1] + strides[0]);
case 20: // { -1, 1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] + strides[1] - strides[0]);
case 21: // { 1, 1, -1 }
return SortIndicesPortal.Get(meshIndex - strides[2] + strides[1] + strides[0]);
case 22: // { -1, -1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] - strides[1] - strides[0]);
case 23: // { 1, -1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] - strides[1] + strides[0]);
case 24: // { -1, 1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] + strides[1] - strides[0]);
case 25: // { 1, 1, 1 }
return SortIndicesPortal.Get(meshIndex + strides[2] + strides[1] + strides[0]);
default:
assert(false);
VTKM_ASSERT(false);
// TODO/FIXME: Should probaly return an invalid value or throw an exception instead
return meshIndex; // Need to error out here
}
} // GetNeighbourIndex
@ -251,13 +234,11 @@ public:
// convert to a sort index
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
vtkm::Id slice = this->VertexSlice(meshIndex);
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((slice == 0) ? FrontBit : 0) |
((slice == this->NumSlices - 1) ? BackBit : 0) | ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id3 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0) | ((pos[2] == 0) ? FrontBit : 0) |
((pos[2] == this->MeshSize[2] - 1) ? BackBit : 0);
// in what follows, the boundary conditions always reset wasAscent
// loop downwards so that we pick the same edges as previous versions
@ -289,13 +270,11 @@ public:
// convert to a sort index
vtkm::Id meshIndex = SortOrderPortal.Get(sortIndex);
vtkm::Id slice = this->VertexSlice(meshIndex);
vtkm::Id row = this->VertexRow(meshIndex);
vtkm::Id col = this->VertexColumn(meshIndex);
vtkm::Int8 boundaryConfig = ((slice == 0) ? FrontBit : 0) |
((slice == this->NumSlices - 1) ? BackBit : 0) | ((col == 0) ? LeftBit : 0) |
((col == this->NumColumns - 1) ? RightBit : 0) | ((row == 0) ? TopBit : 0) |
((row == this->NumRows - 1) ? BottomBit : 0);
vtkm::Id3 pos = this->VertexPos(meshIndex);
vtkm::Int8 boundaryConfig = ((pos[0] == 0) ? LeftBit : 0) |
((pos[0] == this->MeshSize[0] - 1) ? RightBit : 0) | ((pos[1] == 0) ? TopBit : 0) |
((pos[1] == this->MeshSize[1] - 1) ? BottomBit : 0) | ((pos[2] == 0) ? FrontBit : 0) |
((pos[2] == this->MeshSize[2] - 1) ? BackBit : 0);
// Initialize "union find"
int parentId[N_ALL_NEIGHBOURS];

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ArcComparator.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ArcComparator.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CMakeLists.txt → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CMakeLists.txt
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedOtherStartIndexNNeighboursWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedOtherStartIndexNNeighboursWorklet.h
View File

2
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedSimulatedSimplicityIndexComparator.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedSimulatedSimplicityIndexComparator.h
View File

@ -65,7 +65,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVector.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVector.h>
namespace vtkm
{

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVector.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVector.h
View File

2
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVectorDifferentFromNext.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVectorDifferentFromNext.h
View File

@ -66,7 +66,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CombinedVector.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVector.h>
namespace vtkm
{

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/CompressNeighboursWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CompressNeighboursWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ComputeMaxNeighboursWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ComputeMaxNeighboursWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/FindStartIndexWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindStartIndexWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/MergeCombinedOtherStartIndexWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeCombinedOtherStartIndexWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/SubtractAssignWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/SubtractAssignWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/contourtreemesh/UpdateCombinedNeighboursWorklet.h → vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/UpdateCombinedNeighboursWorklet.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/freudenthal_2D/Types.h → vtkm/worklet/contourtree_augmented/meshtypes/freudenthal_2D/Types.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/freudenthal_3D/Types.h → vtkm/worklet/contourtree_augmented/meshtypes/freudenthal_3D/Types.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/marchingcubes_3D/Types.h → vtkm/worklet/contourtree_augmented/meshtypes/marchingcubes_3D/Types.h
View File

0
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/CMakeLists.txt → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/CMakeLists.txt
View File

24
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundary2D.h → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundary2D.h
View File

@ -88,34 +88,34 @@ public:
const InFieldPortalType sortIndicesPortal,
const MeshBoundaryType& meshBoundary,
vtkm::Id& boundaryVertex,
vtkm::Id& boundarySortIndex)
vtkm::Id& boundarySortIndex) const
{
auto meshStructure2D = meshBoundary.GetMeshStructure();
vtkm::Id numBoundary = 2 * meshStructure2D.NumRows + 2 * meshStructure2D.NumColumns - 4;
vtkm::Id numBoundary = 2 * meshStructure2D.MeshSize[1] + 2 * meshStructure2D.MeshSize[0] - 4;
// For comments: [0] -> column, [1] -> row
// Define the boundaryVertex result
if (boundaryId < meshStructure2D->NumColumns)
if (boundaryId < meshStructure2D.MeshSize[0])
{
boundaryVertex = meshStructure2D.VertexId(0, boundaryId);
boundaryVertex = meshStructure2D.VertexId(vtkm::Id2{ boundaryId, 0 });
}
// then bottom row
else if (boundaryId > numBoundary - meshStructure2D.NumColumns - 1)
else if (boundaryId > numBoundary - meshStructure2D.MeshSize[0] - 1)
{
boundaryVertex = meshStructure2D.VertexId(
meshStructure2D.NumRows - 1, boundaryId + meshStructure2D.NumColumns - numBoundary);
boundaryVertex = meshStructure2D.VertexId(vtkm::Id2{
boundaryId + meshStructure2D.MeshSize[0] - numBoundary, meshStructure2D.MeshSize[1] - 1 });
}
// then the row ends
else
{ // row ends
vtkm::Id row = ((boundaryId - meshStructure2D.NumColumns) / 2) + 1;
vtkm::Id col =
((boundaryId - meshStructure2D.NumColumns) % 2) ? (meshStructure2D.NumColumnsn - 1) : 0;
boundaryVertex = meshStructure2D.VertexId(row, col);
boundaryVertex = meshStructure2D.VertexId(vtkm::Id2{
((boundaryId - meshStructure2D.MeshSize[0]) % 2) ? (meshStructure2D.MeshSize[0] - 1) : 0,
((boundaryId - meshStructure2D.MeshSize[0]) / 2) + 1 });
} // row ends
// and fill in the sort index array as well
boundarySortIndex = sortIndicesPortal.Get(boundaryVertex);
/*
/* TODO/FIXME: Delete this comment after code review and tests
// compute how many elements are needed
indexType nBoundary = 2 * nRows + 2 * nCols - 4;

55
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundary3D.h → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundary3D.h
View File

@ -88,33 +88,32 @@ public:
const InFieldPortalType sortIndicesPortal,
const MeshBoundaryType& meshBoundary,
vtkm::Id& boundaryVertex,
vtkm::Id& boundarySortIndex)
vtkm::Id& boundarySortIndex) const
{
auto meshStructure3D = meshBoundary.GetMeshStructure();
vtkm::Id nRows = meshStructure3D.NumRows;
vtkm::Id nCols = meshStructure3D.NumCols;
vtkm::Id nSlices = meshStructure3D.NumSlices;
vtkm::Id3 meshSize = meshStructure3D.MeshSize;
// for comments [0]/x -> column, [1]/y -> row, [2]/z -> slice
// calculate the number of boundary elements - all of the two xy faces
vtkm::Id nBoundary = 2 * nRows * nCols // xy faces
+ 2 * nRows * (nSlices - 2) // yz faces - excluding vertices on xy
+ 2 * (nCols - 2) * (nSlices - 2); // xz face interiors
vtkm::Id nBoundary = 2 * meshSize[1] * meshSize[0] // xy faces
+ 2 * meshSize[1] * (meshSize[2] - 2) // yz faces - excluding vertices on xy
+ 2 * (meshSize[0] - 2) * (meshSize[2] - 2); // xz face interiors
vtkm::Id row = 0, col = 0, slice = 0;
vtkm::Id sliceSize = nRows * nCols;
vtkm::Id sliceBoundarySize = 2 * nRows + 2 * nCols - 4;
vtkm::Id3 pos{ 0, 0, 0 };
vtkm::Id sliceSize = meshSize[1] * meshSize[0];
vtkm::Id sliceBoundarySize = 2 * meshSize[1] + 2 * meshSize[0] - 4;
// do top plane first
if (boundaryId < sliceSize)
{ // top plane
row = boundaryId / nCols;
col = boundaryId % nCols;
slice = 0;
pos[0] = boundaryId % meshSize[0];
pos[1] = boundaryId / meshSize[0];
pos[2] = 0;
} // top plane
// then bottom plane
else if (boundaryId >= nBoundary - sliceSize)
{ // bottom plane
row = (boundaryId - (nBoundary - sliceSize)) / nCols;
col = (boundaryId - (nBoundary - sliceSize)) % nCols;
slice = nSlices - 1;
pos[0] = (boundaryId - (nBoundary - sliceSize)) % meshSize[0];
pos[1] = (boundaryId - (nBoundary - sliceSize)) / meshSize[0];
pos[2] = meshSize[2] - 1;
} // bottom plane
// now we have to deal with the exterior of the remaining slices
else
@ -122,28 +121,28 @@ public:
// first we subtract the size of the first slice
vtkm::Id offsetBoundaryid = boundaryId - sliceSize;
// now we can compute the slice id
slice = 1 + offsetBoundaryid / sliceBoundarySize;
pos[2] = 1 + offsetBoundaryid / sliceBoundarySize;
// compute the local id on the slice
vtkm::Id sliceBoundaryid = offsetBoundaryid % sliceBoundarySize;
// now test for the first and last row
if (sliceBoundaryid < nCols)
vtkm::Id sliceBoundaryId = offsetBoundaryid % sliceBoundarySize;
// now test for the first and last pos[1]
if (sliceBoundaryId < meshSize[0])
{ // first row
row = 0;
col = sliceBoundaryid;
pos[0] = sliceBoundaryId;
pos[1] = 0;
} // first row
else if (sliceBoundaryid >= (sliceBoundarySize - nCols))
else if (sliceBoundaryId >= (sliceBoundarySize - meshSize[0]))
{ // last row
row = nRows - 1;
col = sliceBoundaryid - (sliceBoundarySize - nCols);
pos[0] = sliceBoundaryId - (sliceBoundarySize - meshSize[0]);
pos[1] = meshSize[1] - 1;
} // last row
else
{ // any other row
row = ((sliceBoundaryid - nCols) / 2) + 1;
col = ((sliceBoundaryid - nCols) % 2) ? (nCols - 1) : 0;
pos[0] = ((sliceBoundaryId - meshSize[0]) % 2) ? (meshSize[0] - 1) : 0;
pos[1] = ((sliceBoundaryId - meshSize[0]) / 2) + 1;
} // any other row
} // slice exteriors
// now we have row, col, slice all set, compute the actual ID
boundaryVertex = meshStructure3D.VertexId(slice, row, col);
boundaryVertex = meshStructure3D.VertexId(pos);
// and fill in the index array as well
boundarySortIndex = sortIndicesPortal.Get(boundaryVertex);

6
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h
View File

@ -81,12 +81,12 @@ public:
VTKM_EXEC_CONT
ComputeMeshBoundaryContourTreeMesh() {}
template <typename InFieldPortalType, typename MeshBoundaryType>
template <typename MeshBoundaryType>
VTKM_EXEC void operator()(const vtkm::Id& nodeIndex,
const MeshBoundaryType& meshBoundary,
vtkm::Id& isOnBoundary)
bool& isOnBoundary) const
{
isOnBoundary = meshBoundary.liesOnBoundary(nodeIndex);
isOnBoundary = meshBoundary.LiesOnBoundary(nodeIndex);
/*
indexVector isOnBoundary(globalMeshIndex.size());
for (indexType node = 0; node < globalMeshIndex.size(); node++)

94
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundary2D.h → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary2D.h
View File

@ -62,7 +62,7 @@
#include <cstdlib>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/MeshStructure2D.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure2D.h>
#include <vtkm/cont/ExecutionObjectBase.h>
@ -79,55 +79,96 @@ class MeshBoundary2D
{
public:
// Sort indicies types
using SortOrderPortalType = typename IdArrayType::template ExecutionTypes<DeviceTag>::PortalConst;
using SortIndicesPortalType =
typename IdArrayType::template ExecutionTypes<DeviceTag>::PortalConst;
VTKM_EXEC_CONT
MeshBoundary2D()
: MeshStructure(mesh_dem::MeshStructure2D<DeviceTag>(0, 0))
: MeshStructure(0, 0)
{
}
VTKM_CONT
MeshBoundary2D(vtkm::Id nrows,
vtkm::Id ncols,
const IdArrayType& sortOrder,
vtkm::cont::Token& token)
: MeshStructure(mesh_dem::MeshStructure2D<DeviceTag>(nrows, ncols))
MeshBoundary2D(vtkm::Id2 meshSize, const IdArrayType& sortIndices, vtkm::cont::Token& token)
: MeshStructure(meshSize)
{
this->SortOrderPortal = sortOrder.PrepareForInput(DeviceTag(), token);
this->SortIndicesPortal = sortIndices.PrepareForInput(DeviceTag(), token);
}
VTKM_EXEC_CONT
bool liesOnBoundary(const vtkm::Id index) const
bool LiesOnBoundary(const vtkm::Id meshIndex) const
{
vtkm::Id meshSortOrderValue = this->SortOrderPortal.Get(index);
const vtkm::Id row = this->MeshStructure.VertexRow(meshSortOrderValue);
const vtkm::Id col = this->MeshStructure.VertexColumn(meshSortOrderValue);
return (row == 0) || (col == 0) || (row == this->MeshStructure.NumRows - 1) ||
(col == this->MeshStructure.NumColumns - 1);
const vtkm::Id2 pos{ this->MeshStructure.VertexPos(meshIndex) };
return (pos[0] == 0) || (pos[1] == 0) || (pos[0] == this->MeshStructure.MeshSize[0] - 1) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1);
}
VTKM_EXEC_CONT
const mesh_dem::MeshStructure2D<DeviceTag>& GetMeshStructure() const
bool IsNecessary(const vtkm::Id meshIndex) const
{
vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
const vtkm::Id2 pos{ this->MeshStructure.VertexPos(meshIndex) };
// Keep only when lying on boundary
if ((pos[0] == 0) || (pos[1] == 0) || (pos[0] == this->MeshStructure.MeshSize[0] - 1) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1))
{
// If vertex lies on boundary, keep all corners in any case
if (((pos[1] == 0) && ((pos[0] == 0) || (pos[0] == this->MeshStructure.MeshSize[0] - 1))) ||
((pos[1] == this->MeshStructure.MeshSize[1] - 1) &&
((pos[0] == 0) || (pos[0] == this->MeshStructure.MeshSize[0] - 1))))
{
return true;
}
else
{
// if not a corner, keeep only vertices that are local extrema
vtkm::Id sp, sn;
if (pos[1] == 0 || pos[1] == this->MeshStructure.MeshSize[1] - 1)
{
assert(pos[0] > 0 && pos[0] < this->MeshStructure.MeshSize[0] - 1);
assert(meshIndex >= 1);
sp = this->SortIndicesPortal.Get(meshIndex - 1);
assert(meshIndex + 1 < this->SortIndicesPortal.GetNumberOfValues());
sn = this->SortIndicesPortal.Get(meshIndex + 1);
}
else if (pos[0] == 0 || pos[0] == this->MeshStructure.MeshSize[0] - 1)
{
assert(pos[1] > 0 && pos[1] < this->MeshStructure.MeshSize[1] - 1);
assert(meshIndex >= this->MeshStructure.MeshSize[0]);
sp = this->SortIndicesPortal.Get(meshIndex - this->MeshStructure.MeshSize[0]);
assert(meshIndex + this->MeshStructure.MeshSize[0] <
this->SortIndices.GetNumberOfValues());
sn = this->SortIndicesPortal.Get(meshIndex + this->MeshStructure.MeshSize[0]);
}
return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
}
}
else
{
// Discard vertices in the interior
return false;
}
}
VTKM_EXEC_CONT
const data_set_mesh::MeshStructure2D<DeviceTag>& GetMeshStructure() const
{
return this->MeshStructure;
}
private:
// 2D Mesh size parameters
mesh_dem::MeshStructure2D<DeviceTag> MeshStructure;
SortOrderPortalType SortOrderPortal;
data_set_mesh::MeshStructure2D<DeviceTag> MeshStructure;
SortIndicesPortalType SortIndicesPortal;
};
class MeshBoundary2DExec : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_EXEC_CONT
MeshBoundary2DExec(vtkm::Id nrows, vtkm::Id ncols, const IdArrayType& inSortOrder)
: NumRows(nrows)
, NumColumns(ncols)
, SortOrder(inSortOrder)
MeshBoundary2DExec(vtkm::Id2 inMeshSize, const IdArrayType& inSortIndices)
: MeshSize(inMeshSize)
, SortIndices(inSortIndices)
{
}
@ -135,14 +176,13 @@ public:
template <typename DeviceTag>
MeshBoundary2D<DeviceTag> PrepareForExecution(DeviceTag, vtkm::cont::Token& token) const
{
return MeshBoundary2D<DeviceTag>(this->NumRows, this->NumColumns, this->SortOrder, token);
return MeshBoundary2D<DeviceTag>(this->MeshSize, this->SortIndices, token);
}
private:
// 2D Mesh size parameters
vtkm::Id NumRows;
vtkm::Id NumColumns;
const IdArrayType& SortOrder;
vtkm::Id2 MeshSize;
const IdArrayType& SortIndices;
};

295
vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary3D.h Normal file
View File

@ -0,0 +1,295 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt 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) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
// This header contains a collection of classes used to describe the boundary
// of a mesh, for each main mesh type (i.e., 2D, 3D, and ContourTreeMesh).
// For each mesh type, there are two classes, the actual boundary desriptor
// class and an ExectionObject class with the PrepareForInput function that
// VTKm expects to generate the object for the execution environment.
#ifndef vtk_m_worklet_contourtree_augmented_mesh_boundary_mesh_boundary_3d_h
#define vtk_m_worklet_contourtree_augmented_mesh_boundary_mesh_boundary_3d_h
#include <cstdlib>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure3D.h>
#include <vtkm/cont/ExecutionObjectBase.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
template <typename DeviceTag>
class MeshBoundary3D : public vtkm::cont::ExecutionObjectBase
{
public:
// Sort indicies types
using SortIndicesPortalType =
typename IdArrayType::template ExecutionTypes<DeviceTag>::PortalConst;
VTKM_EXEC_CONT
MeshBoundary3D()
: MeshStructure(data_set_mesh::MeshStructure3D<DeviceTag>(vtkm::Id3{ 0, 0, 0 }))
{
}
VTKM_CONT
MeshBoundary3D(vtkm::Id3 meshSize, const IdArrayType& inSortIndices, vtkm::cont::Token& token)
: MeshStructure(data_set_mesh::MeshStructure3D<DeviceTag>(meshSize))
{
this->SortIndicesPortal = inSortIndices.PrepareForInput(DeviceTag(), token);
}
VTKM_EXEC_CONT
bool LiesOnBoundary(const vtkm::Id meshIndex) const
{
const vtkm::Id3 pos = this->MeshStructure.VertexPos(meshIndex);
return (pos[0] == 0) || (pos[1] == 0) || (pos[2] == 0) ||
(pos[0] == this->MeshStructure.MeshSize[0] - 1) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1) ||
(pos[2] == this->MeshStructure.MeshSize[2] - 1);
}
VTKM_EXEC_CONT
vtkm::Id CountLinkComponentsIn2DSlice(const vtkm::Id meshIndex, const vtkm::Id2 strides) const
{
// IMPORTANT: We assume that function is called only for *interior* vertices (i.e., neither row nor col
// within slice is 0 and we do not need to check for boundary cases).
vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
bool prevWasInUpperLink = false;
vtkm::Id numComponents = 0;
const int N_INCIDENT_EDGES_2D = 6;
for (vtkm::Id edgeNo = 0; edgeNo < N_INCIDENT_EDGES_2D; edgeNo++)
{ // per edge
VTKM_ASSERT(meshIndex + strides[1] + strides[0] <
this->SortIndicesPortal.GetNumberOfValues());
VTKM_ASSERT(meshIndex - strides[1] - strides[0] >= 0);
vtkm::Id nbrSortIndex;
switch (edgeNo)
{
case 0:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[0]);
break; // [1] , [0] + 1
case 1:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1] + strides[0]);
break; // [1] + 1, [0] + 1
case 2:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1]);
break; // [1] + 1, [0]
case 3:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[0]);
break; // [1] , [0] - 1
case 4:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1] - strides[0]);
break; // [1] - 1, [0] - 1
case 5:
nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1]);
break; // [1] - 1, [0]
default:
std::abort();
}
bool currIsInUpperLink = (nbrSortIndex > sortIndex);
numComponents += (edgeNo != 0 && currIsInUpperLink != prevWasInUpperLink) ? 1 : 0;
} // per edge
return numComponents;
}
VTKM_EXEC_CONT
bool IsNecessary(vtkm::Id meshIndex) const
{
vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
vtkm::Id3 pos{ this->MeshStructure.VertexPos(meshIndex) };
vtkm::Id nPerSlice = this->MeshStructure.MeshSize[0] *
this->MeshStructure.MeshSize[1]; // number of vertices on a [2]-perpendicular "slice"
// Keep only when lying on boundary
if ((pos[1] == 0) || (pos[0] == 0) || (pos[2] == 0) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1) ||
(pos[0] == this->MeshStructure.MeshSize[0] - 1) ||
(pos[2] == this->MeshStructure.MeshSize[2] - 1))
{
// Keep data on corners
bool atEndOfLine = (pos[0] == 0) || (pos[0] == this->MeshStructure.MeshSize[0] - 1);
bool atQuadCorner = (pos[1] == 0 && atEndOfLine) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1 && atEndOfLine);
if ((pos[2] == 0 && atQuadCorner) ||
(pos[2] == this->MeshStructure.MeshSize[2] - 1 && atQuadCorner))
{
return true;
}
else
{
// Check if vertex lies along one of the boundary edges, if so, keep
// local extrema
// Edges in [0] direction
if ((pos[1] == 0 && pos[2] == 0) ||
(pos[1] == 0 && pos[2] == this->MeshStructure.MeshSize[2] - 1) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1 && pos[2] == 0) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1 &&
pos[2] == this->MeshStructure.MeshSize[2] - 1))
{
VTKM_ASSERT(meshIndex >= 1);
vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - 1);
VTKM_ASSERT(meshIndex + 1 < this->SortIndicesPortal.GetNumberOfValues());
vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + 1);
return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
}
// Edges in [1] directtion
else if ((pos[0] == 0 && pos[2] == 0) ||
(pos[0] == 0 && pos[2] == this->MeshStructure.MeshSize[2] - 1) ||
(pos[0] == this->MeshStructure.MeshSize[0] - 1 && pos[2] == 0) ||
(pos[0] == this->MeshStructure.MeshSize[0] - 1 &&
pos[2] == this->MeshStructure.MeshSize[2] - 1))
{
VTKM_ASSERT(pos[1] > 0 && pos[1] < this->MeshStructure.MeshSize[1] - 1);
VTKM_ASSERT(meshIndex >= this->MeshStructure.MeshSize[0]);
vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - this->MeshStructure.MeshSize[0]);
VTKM_ASSERT(meshIndex + this->MeshStructure.MeshSize[0] <
this->SortIndicesPortal.GetNumberOfValues());
vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + this->MeshStructure.MeshSize[0]);
return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
}
// Edges in [2] direction
else if ((pos[1] == 0 && pos[0] == 0) ||
(pos[1] == 0 && pos[0] == this->MeshStructure.MeshSize[0] - 1) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1 && pos[0] == 0) ||
(pos[1] == this->MeshStructure.MeshSize[1] - 1 &&
pos[0] == this->MeshStructure.MeshSize[0] - 1))
{
VTKM_ASSERT(meshIndex >= nPerSlice);
vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - nPerSlice);
VTKM_ASSERT(meshIndex + nPerSlice < this->SortIndicesPortal.GetNumberOfValues());
vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + nPerSlice);
return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
}
else
{
// On a face/slice
if (pos[2] == 0 || pos[2] == this->MeshStructure.MeshSize[2] - 1)
{ // On [2]-perpendicular face
VTKM_ASSERT(pos[0] != 0 && pos[0] != this->MeshStructure.MeshSize[0]);
VTKM_ASSERT(pos[1] != 0 && pos[1] != this->MeshStructure.MeshSize[1]);
return CountLinkComponentsIn2DSlice(meshIndex, this->MeshStructure.MeshSize[0], 1) ==
1; // FIXME: or != 2;
}
else if (pos[1] == 0 || pos[1] == this->MeshStructure.MeshSize[1] - 1)
{ // On [1]-perpendicular face
VTKM_ASSERT(pos[0] != 0 && pos[0] != this->MeshStructure.MeshSize[0]);
VTKM_ASSERT(pos[2] != 0 && pos[2] != this->MeshStructure.MeshSize[2]);
return CountLinkComponentsIn2DSlice(meshIndex, nPerSlice, 1) == 1; // FIXME: or != 2;
}
else
{ // On [0]-perpendicular face
VTKM_ASSERT(pos[0] == 0 || pos[0] == this->MeshStructure.MeshSize[0] - 1);
VTKM_ASSERT(pos[1] != 0 && pos[1] != this->MeshStructure.MeshSize[1]);
VTKM_ASSERT(pos[2] != 0 && pos[2] != this->MeshStructure.MeshSize[2]);
return CountLinkComponentsIn2DSlice(
meshIndex, this->MeshStructure.MeshSize[0], nPerSlice) == 1; // FIXME: or != 2;
}
}
}
}
else
{
return false;
}
}
VTKM_EXEC_CONT
const data_set_mesh::MeshStructure3D<DeviceTag>& GetMeshStructure() const
{
return this->MeshStructure;
}
protected:
// 3D Mesh size parameters
data_set_mesh::MeshStructure3D<DeviceTag> MeshStructure;
SortIndicesPortalType SortIndicesPortal;
};
class MeshBoundary3DExec : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_EXEC_CONT
MeshBoundary3DExec(vtkm::Id3 meshSize, const IdArrayType& inSortIndices)
: MeshSize(meshSize)
, SortIndices(inSortIndices)
{
}
VTKM_CONT
template <typename DeviceTag>
MeshBoundary3D<DeviceTag> PrepareForExecution(DeviceTag, vtkm::cont::Token& token) const
{
return MeshBoundary3D<DeviceTag>(this->MeshSize, this->SortIndices, token);
}
protected:
// 3D Mesh size parameters
vtkm::Id3 MeshSize;
const IdArrayType& SortIndices;
};
} // namespace contourtree_augmented
} // worklet
} // vtkm
#endif

46
vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h → vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h
View File

@ -84,32 +84,35 @@ public:
VTKM_CONT
MeshBoundaryContourTreeMesh(const IdArrayType& globalMeshIndex,
vtkm::Id totalNRows,
vtkm::Id totalNCols,
vtkm::Id3 globalSize,
vtkm::Id3 minIdx,
vtkm::Id3 maxIdx,
vtkm::cont::Token& token)
: TotalNRows(totalNRows)
, TotalNCols(totalNCols)
: GlobalSize(globalSize)
, MinIdx(minIdx)
, MaxIdx(maxIdx)
{
assert(this->TotalNRows > 0 && this->TotalNCols > 0);
assert(this->GlobalSize[0] > 0 && this->GlobalSize[1] > 0);
this->GlobalMeshIndexPortal = globalMeshIndex.PrepareForInput(DeviceTag(), token);
}
VTKM_EXEC_CONT
bool liesOnBoundary(const vtkm::Id index) const
bool LiesOnBoundary(const vtkm::Id index) const
{
vtkm::Id idx = this->GlobalMeshIndexPortal.Get(index);
vtkm::Id3 rcs;
rcs[0] = vtkm::Id((idx % (this->TotalNRows * this->TotalNCols)) / this->TotalNCols);
rcs[1] = vtkm::Id(idx % this->TotalNCols);
rcs[2] = vtkm::Id(idx / (this->TotalNRows * this->TotalNCols));
for (int d = 0; d < 3; ++d)
vtkm::Id global_idx = this->GlobalMeshIndexPortal.Get(index);
vtkm::Id3 mesh_idx{ vtkm::Id(global_idx % this->GlobalSize[0]),
vtkm::Id((global_idx % (this->GlobalSize[0] * this->GlobalSize[1])) /
this->GlobalSize[0]),
vtkm::Id(global_idx / (this->GlobalSize[0] * this->GlobalSize[1])) };
// FIXME: Probably better communicate n_dims in constructor or make it a template parameter
// Or at least be more consistent in setting this in MinIdx/MaxIdx. Currently MinIdx[2] is 0
// and MaxIdx[2] is -1 for a 2D data set.
const auto n_dims = (MaxIdx[2] == -1) ? 2 : 3;
for (int d = 0; d < n_dims; ++d)
{
if (this->MinIdx[d] != this->MaxIdx[d] &&
(rcs[d] == this->MinIdx[d] || rcs[d] == this->MaxIdx[d]))
(mesh_idx[d] == this->MinIdx[d] || mesh_idx[d] == this->MaxIdx[d]))
{
return true;
}
@ -117,10 +120,12 @@ public:
return false;
}
VTKM_EXEC_CONT
bool IsNecessary(const vtkm::Id idx) const { return this->LiesOnBoundary(idx); }
private:
// mesh block parameters
vtkm::Id TotalNRows;
vtkm::Id TotalNCols;
vtkm::Id3 GlobalSize;
vtkm::Id3 MinIdx;
vtkm::Id3 MaxIdx;
IndicesPortalType GlobalMeshIndexPortal;
@ -132,13 +137,11 @@ class MeshBoundaryContourTreeMeshExec : public vtkm::cont::ExecutionObjectBase
public:
VTKM_EXEC_CONT
MeshBoundaryContourTreeMeshExec(const IdArrayType& globalMeshIndex,
vtkm::Id totalNRows,
vtkm::Id totalNCols,
vtkm::Id3 globalSize,
vtkm::Id3 minIdx,
vtkm::Id3 maxIdx)
: GlobalMeshIndex(globalMeshIndex)
, TotalNRows(totalNRows)
, TotalNCols(totalNCols)
, GlobalSize(globalSize)
, MinIdx(minIdx)
, MaxIdx(maxIdx)
{
@ -150,13 +153,12 @@ public:
vtkm::cont::Token& token) const
{
return MeshBoundaryContourTreeMesh<DeviceTag>(
this->GlobalMeshIndex, this->TotalNRows, this->TotalNCols, this->MinIdx, this->MaxIdx, token);
this->GlobalMeshIndex, this->GlobalSize, this->MinIdx, this->MaxIdx, token);
}
private:
const IdArrayType& GlobalMeshIndex;
vtkm::Id TotalNRows;
vtkm::Id TotalNCols;
vtkm::Id3 GlobalSize;
vtkm::Id3 MinIdx;
vtkm::Id3 MaxIdx;
};

15
vtkm/worklet/contourtree_distributed/MergeBlockFunctor.h
View File

@ -54,6 +54,7 @@
#define vtk_m_worklet_contourtree_distributed_mergeblockfunctor_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
// clang-format off
@ -117,7 +118,8 @@ void MergeBlockFunctor(
// Construct the two contour tree mesh by assignign the block data
vtkm::worklet::contourtree_augmented::ContourTreeMesh<FieldType> contourTreeMeshIn;
contourTreeMeshIn.NumVertices = recvblock.NumVertices;
contourTreeMeshIn.SortOrder = recvblock.SortOrder;
contourTreeMeshIn.SortOrder = vtkm::cont::ArrayHandleIndex(contourTreeMeshIn.NumVertices);
contourTreeMeshIn.SortIndices = vtkm::cont::ArrayHandleIndex(contourTreeMeshIn.NumVertices);
contourTreeMeshIn.SortedValues = recvblock.SortedValue;
contourTreeMeshIn.GlobalMeshIndex = recvblock.GlobalMeshIndex;
contourTreeMeshIn.Neighbours = recvblock.Neighbours;
@ -126,7 +128,8 @@ void MergeBlockFunctor(
vtkm::worklet::contourtree_augmented::ContourTreeMesh<FieldType> contourTreeMeshOut;
contourTreeMeshOut.NumVertices = block->NumVertices;
contourTreeMeshOut.SortOrder = block->SortOrder;
contourTreeMeshOut.SortOrder = vtkm::cont::ArrayHandleIndex(contourTreeMeshOut.NumVertices);
contourTreeMeshOut.SortIndices = vtkm::cont::ArrayHandleIndex(contourTreeMeshOut.NumVertices);
contourTreeMeshOut.SortedValues = block->SortedValue;
contourTreeMeshOut.GlobalMeshIndex = block->GlobalMeshIndex;
contourTreeMeshOut.Neighbours = block->Neighbours;
@ -158,7 +161,6 @@ void MergeBlockFunctor(
{
// Save the data from our block for the next iteration
block->NumVertices = contourTreeMeshOut.NumVertices;
block->SortOrder = contourTreeMeshOut.SortOrder;
block->SortedValue = contourTreeMeshOut.SortedValues;
block->GlobalMeshIndex = contourTreeMeshOut.GlobalMeshIndex;
block->Neighbours = contourTreeMeshOut.Neighbours;
@ -180,11 +182,7 @@ void MergeBlockFunctor(
currBlockOrigin[1] + currBlockSize[1] - 1,
currBlockOrigin[2] + currBlockSize[2] - 1);
auto meshBoundaryExecObj =
contourTreeMeshOut.GetMeshBoundaryExecutionObject(globalSize[0], // totalNRows
globalSize[1], // totalNCols
currBlockOrigin, // minIdx
maxIdx // maxIdx
);
contourTreeMeshOut.GetMeshBoundaryExecutionObject(globalSize, currBlockOrigin, maxIdx);
worklet.Run(
contourTreeMeshOut.SortedValues, // Unused param. Provide something to keep the API happy
contourTreeMeshOut,
@ -215,7 +213,6 @@ void MergeBlockFunctor(
// Copy the data from newContourTreeMesh into block
block->NumVertices = newContourTreeMesh->NumVertices;
block->SortOrder = newContourTreeMesh->SortOrder;
block->SortedValue = newContourTreeMesh->SortedValues;
block->GlobalMeshIndex = newContourTreeMesh->GlobalMeshIndex;
block->Neighbours = newContourTreeMesh->Neighbours;

15
vtkm/worklet/contourtree_distributed/MultiBlockContourTreeHelper.h
View File

@ -56,7 +56,7 @@
#include <vtkm/worklet/contourtree_distributed/SpatialDecomposition.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem_meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/Types.h>
#include <vtkm/cont/BoundsCompute.h>
@ -64,7 +64,7 @@
//#include <vtkm/cont/AssignerPartitionedDataSet.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/worklet/contourtree_augmented/mesh_dem/IdRelabeler.h>
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h>
namespace vtkm
@ -153,13 +153,6 @@ public:
unsigned int computeRegularStructure)
{
vtkm::Id startRow = localBlockOrigin[0];
vtkm::Id startCol = localBlockOrigin[1];
vtkm::Id startSlice = localBlockOrigin[2];
vtkm::Id numRows = localBlockSize[0];
vtkm::Id numCols = localBlockSize[1];
vtkm::Id totalNumRows = globalSize[0];
vtkm::Id totalNumCols = globalSize[1];
// compute the global mesh index and initalize the local contour tree mesh
if (computeRegularStructure == 1)
{
@ -170,7 +163,7 @@ public:
vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler>(
sortOrder,
vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler(
startRow, startCol, startSlice, numRows, numCols, totalNumRows, totalNumCols));
localBlockOrigin, localBlockSize, globalSize));
vtkm::cont::Algorithm::Copy(transformedIndex, localGlobalMeshIndex);
// Compute the local contour tree mesh
auto localContourTreeMesh = new vtkm::worklet::contourtree_augmented::ContourTreeMesh<T>(
@ -190,7 +183,7 @@ public:
auto transformedIndex = vtkm::cont::make_ArrayHandleTransform(
permutedSortOrder,
vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler(
startRow, startCol, startSlice, numRows, numCols, totalNumRows, totalNumCols));
localBlockOrigin, localBlockSize, globalSize));
vtkm::cont::Algorithm::Copy(transformedIndex, localGlobalMeshIndex);
// Compute the local contour tree mesh
auto localContourTreeMesh = new vtkm::worklet::contourtree_augmented::ContourTreeMesh<T>(

41
vtkm/worklet/testing/UnitTestContourTreeUniformAugmented.cxx
View File

@ -310,19 +310,17 @@ private:
template <typename FieldType, typename StorageType>
void CallTestContourTreeAugmentedSteps(
const vtkm::cont::ArrayHandle<FieldType, StorageType> fieldArray,
const vtkm::Id nRows,
const vtkm::Id nCols,
const vtkm::Id nSlices,
const vtkm::Id3 meshSize,
bool useMarchingCubes,
unsigned int computeRegularStructure,
ExpectedStepResults& expectedResults) const
{
using namespace vtkm::worklet::contourtree_augmented;
// 2D Contour Tree
if (nSlices == 1)
if (meshSize[2] == 1)
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_2D_Freudenthal<FieldType, StorageType> mesh(nRows, nCols);
DataSetMeshTriangulation2DFreudenthal mesh(vtkm::Id2{ meshSize[0], meshSize[1] });
// Run the contour tree on the mesh
RunTestContourTreeAugmentedSteps(fieldArray,
mesh,
@ -335,7 +333,7 @@ private:
else if (useMarchingCubes)
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_3D_MarchingCubes<FieldType, StorageType> mesh(nRows, nCols, nSlices);
DataSetMeshTriangulation3DMarchingCubes mesh(meshSize);
// Run the contour tree on the mesh
RunTestContourTreeAugmentedSteps(fieldArray,
mesh,
@ -348,7 +346,7 @@ private:
else
{
// Build the mesh and fill in the values
Mesh_DEM_Triangulation_3D_Freudenthal<FieldType, StorageType> mesh(nRows, nCols, nSlices);
DataSetMeshTriangulation3DFreudenthal mesh(meshSize);
// Run the contour tree on the mesh
RunTestContourTreeAugmentedSteps(fieldArray,
mesh,
@ -377,16 +375,13 @@ private:
dataSet.GetCellSet().CopyTo(cellSet);
vtkm::Id3 pointDimensions = cellSet.GetPointDimensions();
vtkm::Id nRows = pointDimensions[0];
vtkm::Id nCols = pointDimensions[1];
vtkm::Id nSlices = pointDimensions[2];
vtkm::cont::ArrayHandle<vtkm::Float32> field;
dataSet.GetField("pointvar").GetData().CopyTo(field);
// Run the specific test
CallTestContourTreeAugmentedSteps(
field, nRows, nCols, nSlices, useMarchingCubes, computeRegularStructure, expectedResults);
field, pointDimensions, useMarchingCubes, computeRegularStructure, expectedResults);
}
///
@ -727,10 +722,8 @@ public:
vtkm::cont::CellSetStructured<2> cellSet;
dataSet.GetCellSet().CopyTo(cellSet);
vtkm::Id2 pointDimensions = cellSet.GetPointDimensions();
vtkm::Id nRows = pointDimensions[0];
vtkm::Id nCols = pointDimensions[1];
vtkm::Id nSlices = 1;
vtkm::Id2 pointDimensions2D = cellSet.GetPointDimensions();
vtkm::Id3 meshSize{ pointDimensions2D[0], pointDimensions2D[1], 1 };
vtkm::cont::ArrayHandle<vtkm::Float32> field;
dataSet.GetField("pointvar").GetData().CopyTo(field);
@ -747,9 +740,7 @@ public:
contourTree,
meshSortOrder,
numIterations,
nRows,
nCols,
nSlices,
meshSize,
useMarchingCubes,
computeRegularStructure);
@ -799,9 +790,6 @@ public:
dataSet.GetCellSet().CopyTo(cellSet);
vtkm::Id3 pointDimensions = cellSet.GetPointDimensions();
vtkm::Id nRows = pointDimensions[0];
vtkm::Id nCols = pointDimensions[1];
vtkm::Id nSlices = pointDimensions[2];
vtkm::cont::ArrayHandle<vtkm::Float32> field;
dataSet.GetField("pointvar").GetData().CopyTo(field);
@ -818,9 +806,7 @@ public:
contourTree,
meshSortOrder,
numIterations,
nRows,
nCols,
nSlices,
pointDimensions,
useMarchingCubes,
computeRegularStructure);
@ -877,9 +863,6 @@ public:
dataSet.GetCellSet().CopyTo(cellSet);
vtkm::Id3 pointDimensions = cellSet.GetPointDimensions();
vtkm::Id nRows = pointDimensions[0];
vtkm::Id nCols = pointDimensions[1];
vtkm::Id nSlices = pointDimensions[2];
vtkm::cont::ArrayHandle<vtkm::Float32> field;
dataSet.GetField("pointvar").GetData().CopyTo(field);
@ -896,9 +879,7 @@ public:
contourTree,
meshSortOrder,
numIterations,
nRows,
nCols,
nSlices,
pointDimensions,
useMarchingCubes,
computeRegularStructure);