Skip to content

Commit 778db6b

Browse files
abessiariabessiari
committed
Using GetBranchByVolume Filter
1 parent 756ba13 commit 778db6b

3 files changed

Lines changed: 168 additions & 42 deletions

File tree

‎src/libs/vtkh/filters/ContourTree.cpp‎

Lines changed: 153 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -3,13 +3,15 @@
33
#include <vtkh/filters/Recenter.hpp>
44

55
// vtkm includes
6-
#include <vtkm/cont/DeviceAdapter.h>
7-
#include <vtkm/cont/Storage.h>
8-
#include <vtkm/internal/Configure.h>
6+
#include <vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h>
7+
#include <vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.h>
8+
#include <vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h>
9+
#include <vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h>
910
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h>
1011
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/ProcessContourTree.h>
11-
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/processcontourtree/Branch.h>
12-
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/processcontourtree/PiecewiseLinearFunction.h>
12+
#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
13+
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalContourTree.h>
14+
#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/TreeCompiler.h>
1315

1416
#include <vtkh/filters/GhostStripper.hpp>
1517

@@ -198,14 +200,52 @@ void ContourTree::PostExecute()
198200
Filter::PostExecute();
199201
}
200202

203+
struct TopVolumeBranchSaddleIsoValueFunctor
204+
{
205+
vtkm::Id nSelectedBranches;
206+
std::vector<vtkm::Float64> dbranches;
207+
208+
public:
209+
210+
TopVolumeBranchSaddleIsoValueFunctor(vtkm::Id nSelectedBranches) : nSelectedBranches(nSelectedBranches) {
211+
}
212+
void operator()(const vtkm::cont::ArrayHandle<vtkm::Float32> &arr)
213+
{
214+
auto topVolBranchSaddleIsoValue = arr.ReadPortal();
215+
216+
for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch) {
217+
dbranches.push_back(topVolBranchSaddleIsoValue.Get(branch));
218+
}
219+
}
220+
221+
void operator()(const vtkm::cont::ArrayHandle<vtkm::Float64> &arr)
222+
{
223+
auto topVolBranchSaddleIsoValue = arr.ReadPortal();
224+
225+
for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch) {
226+
dbranches.push_back(topVolBranchSaddleIsoValue.Get(branch));
227+
}
228+
}
229+
template <typename T>
230+
void operator()(const T&)
231+
{
232+
throw vtkm::cont::ErrorBadValue("TopVolumeBranchSaddleIsoValueFunctor Expected Float32 or Float64!");
233+
}
234+
235+
vtkm::Float64 Get(vtkm::Id branch)
236+
{
237+
return dbranches[branch];
238+
}
239+
};
240+
201241
struct AnalyzerFunctor
202242
{
203-
vtkm::filter::scalar_topology::ContourTreeAugmented& filter;
243+
vtkm::filter::scalar_topology::ContourTreeAugmented* filter;
204244
vtkh::ContourTree& contourTree;
205245
bool dataFieldIsSorted;
206246

207247
public:
208-
AnalyzerFunctor(vtkh::ContourTree& contourTree, vtkm::filter::scalar_topology::ContourTreeAugmented& filter): contourTree(contourTree), filter(filter) {
248+
AnalyzerFunctor(vtkh::ContourTree& contourTree, vtkm::filter::scalar_topology::ContourTreeAugmented* filter): contourTree(contourTree), filter(filter) {
209249
}
210250

211251
void SetDataFieldIsSorted(bool dataFieldIsSorted) {
@@ -214,12 +254,12 @@ struct AnalyzerFunctor
214254

215255
void operator()(const vtkm::cont::ArrayHandle<vtkm::Float32> &arr) const
216256
{
217-
contourTree.analysis<vtkm::Float32>(filter, dataFieldIsSorted, arr);
257+
contourTree.analysis<vtkm::Float32>(*filter, dataFieldIsSorted, arr);
218258
}
219259

220260
void operator()(const vtkm::cont::ArrayHandle<vtkm::Float64> &arr) const
221261
{
222-
contourTree.analysis<vtkm::Float64>(filter, dataFieldIsSorted, arr);
262+
contourTree.analysis<vtkm::Float64>(*filter, dataFieldIsSorted, arr);
223263
}
224264

225265
template <typename T>
@@ -343,8 +383,65 @@ template<typename DataValueType> void ContourTree::analysis(vtkm::filter::scalar
343383
}
344384
}
345385

386+
void ContourTree::distributed_analysis(const vtkm::cont::PartitionedDataSet& result, int mpi_rank) {
387+
vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter bd_filter;
388+
auto bd_result = bd_filter.Execute(result);
389+
390+
#if 0
391+
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
392+
{
393+
auto ds = bd_result.GetPartition(ds_no);
394+
std::string branchDecompositionFileName = std::string("BranchDecomposition_Rank_") +
395+
std::to_string(static_cast<int>(mpi_rank)) + std::string("_Block_") +
396+
std::to_string(static_cast<int>(ds_no)) + std::string(".txt");
397+
398+
std::ofstream treeStream(branchDecompositionFileName.c_str());
399+
treeStream
400+
<< vtkm::filter::scalar_topology::HierarchicalVolumetricBranchDecomposer::PrintBranches(ds);
401+
}
402+
#endif
403+
404+
vtkm::filter::scalar_topology::SelectTopVolumeContoursFilter tp_filter;
405+
vtkm::Id numBranches = m_levels;
406+
407+
tp_filter.SetSavedBranches(numBranches);
408+
bd_result = tp_filter.Execute(bd_result);
409+
410+
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
411+
{
412+
auto ds = bd_result.GetPartition(ds_no);
413+
auto topVolBranchGRId = ds.GetField("TopVolumeBranchGlobalRegularIds")
414+
.GetData()
415+
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
416+
.ReadPortal();
417+
418+
auto topVolBranchSaddleEpsilon = ds.GetField("TopVolumeBranchSaddleEpsilon")
419+
.GetData()
420+
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
421+
.ReadPortal();
422+
423+
vtkm::Id nSelectedBranches = topVolBranchGRId.GetNumberOfValues();
424+
425+
TopVolumeBranchSaddleIsoValueFunctor iso_functor(nSelectedBranches);
426+
427+
vtkm::cont::CastAndCall(ds.GetField("TopVolumeBranchSaddleIsoValue").GetData(), iso_functor);
428+
429+
430+
vtkm::Float32 eps = 0.00001;
431+
432+
for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch)
433+
{
434+
m_iso_values[branch] = iso_functor.Get(branch) + (eps * topVolBranchSaddleEpsilon.Get(branch));
435+
}
436+
437+
break;
438+
}
439+
440+
std::sort(m_iso_values.begin(), m_iso_values.end());
441+
}
442+
346443
void ContourTree::DoExecute()
347-
{
444+
{
348445
vtkh::DataSet *old_input = this->m_input;
349446
const int before_num_domains = this->m_input->GetNumberOfDomains();
350447

@@ -415,17 +512,20 @@ void ContourTree::DoExecute()
415512
this->m_input->GetDomain(0, dom, domain_id);
416513
inDataSet.AppendPartition(dom);
417514

515+
#ifdef DEBUG
418516
if( mpi_size != 1 )
419517
{
420518
std::ostringstream ostr;
421519
ostr << "rank: " << mpi_rank
422520
<< " coord system range: " << dom.GetCoordinateSystem(0).GetRange() << std::endl;
423521
std::cout << ostr.str();
424522
}
425-
#ifdef DEBUG
426523
#endif
427524
#endif // VTKH_PARALLEL
428525

526+
vtkm::filter::Filter *filterField;
527+
528+
#ifndef VTKH_PARALLEL
429529
bool useMarchingCubes = false;
430530
// Compute the fully augmented contour tree.
431531
// This should always be true for now in order for the isovalue selection to work.
@@ -434,43 +534,60 @@ void ContourTree::DoExecute()
434534
//Convert the mesh of values into contour tree, pairs of vertex ids
435535
vtkm::filter::scalar_topology::ContourTreeAugmented filter(useMarchingCubes, computeRegularStructure);
436536

437-
filter.SetActiveField(m_field_name);
438-
439-
#ifdef VTKH_PARALLEL
537+
#ifdef DEBUG
538+
std::cout << "--- BEGIN_SUMMARY inDataSet" << std::endl;
539+
inDataSet.PrintSummary( std::cout );
540+
std::cout << "--- END_SUMMARY inDataSet" << std::endl;
541+
#endif
542+
#else // VTKH_PARALLEL
543+
// TODO SET VTKM to do vtkm::cont::LogLovel::Info
544+
// vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter(vtkm::cont::LogLevel::Info, vtkm::cont::LogLevel::Info);
545+
vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter;
546+
vtkm::filter::scalar_topology::ContourTreeAugmented aug_filter(false, true);
547+
548+
if (mpi_size == 1) {
549+
filterField = &aug_filter;
550+
} else {
551+
filter.SetUseBoundaryExtremaOnly(true);
552+
filter.SetUseMarchingCubes(false);
553+
filter.SetAugmentHierarchicalTree(true);
440554
ShiftLogicalOriginToZero(inDataSet);
441555
ComputeGlobalPointSize(inDataSet);
556+
filterField = &filter;
557+
}
558+
559+
#ifdef DEBUG
560+
std::cout << "--- BEGIN_SUMMARY inDataSet:mpi_rank=" << mpi_rank << std::endl;
561+
inDataSet.PrintSummary( std::cout );
562+
std::cout << "--- END_SUMMARY inDataSet:mpi_rank=" << mpi_rank << std::endl;
563+
#endif
442564
#endif // VTKH_PARALLEL
443565

444-
auto result = filter.Execute(inDataSet);
566+
filterField->SetActiveField(m_field_name);
445567

446-
m_iso_values.resize(m_levels);
568+
auto result = filterField->Execute(inDataSet);
447569

448-
if (mpi_rank == 0) {
449-
AnalyzerFunctor analyzerFunctor(*this, filter);
570+
m_iso_values.resize(m_levels);
450571

451572
#ifndef VTKH_PARALLEL
573+
if (mpi_rank == 0) {
574+
AnalyzerFunctor analyzerFunctor(*this, &filter);
452575
analyzerFunctor.SetDataFieldIsSorted(false);
453576
vtkm::cont::CastAndCall(inDataSet.GetField(m_field_name).GetData(), analyzerFunctor);
577+
} // mpi_rank == 0
454578
#else
455-
if(mpi_size == 1)
456-
{
457-
analyzerFunctor.SetDataFieldIsSorted(false);
458-
vtkm::cont::CastAndCall(inDataSet.GetPartitions()[0].GetField(m_field_name).GetData(), analyzerFunctor);
459-
} else {
460-
analyzerFunctor.SetDataFieldIsSorted(true);
461-
462-
/*
463-
if( result.GetPartitions()[0].GetNumberOfFields() > 1 ) {
464-
vtkm::cont::CastAndCall(result.GetPartitions()[0].GetField("values").GetData(), analyzerFunctor);
465-
} else {
466-
vtkm::cont::CastAndCall(result.GetPartitions()[0].GetField(0).GetData(), analyzerFunctor);
467-
}*/
468-
469-
// TODO TO BE REVISITED. Tested with: srun -n 8 ./t_vtk-h_contour_tree_par
470-
vtkm::cont::CastAndCall(result.GetPartitions()[0].GetField("resultData").GetData(), analyzerFunctor);
471-
}
579+
if (mpi_size > 1)
580+
{ // BranchDecomposition
581+
distributed_analysis(result, mpi_rank);
582+
}
583+
else
584+
{
585+
AnalyzerFunctor analyzerFunctor(*this, (vtkm::filter::scalar_topology::ContourTreeAugmented *) filterField);
586+
587+
analyzerFunctor.SetDataFieldIsSorted(false);
588+
vtkm::cont::CastAndCall(inDataSet.GetPartitions()[0].GetField(m_field_name).GetData(), analyzerFunctor);
589+
}
472590
#endif // VTKH_PARALLEL
473-
} // mpi_rank == 0
474591

475592
#ifdef VTKH_PARALLEL
476593
MPI_Bcast(&m_iso_values[0], m_levels, MPI_DOUBLE, 0, mpi_comm);

‎src/libs/vtkh/filters/ContourTree.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ class ContourTree : public Filter
2727
private:
2828
friend class AnalyzerFunctor;
2929
template<typename DataValueType> void analysis(vtkm::filter::scalar_topology::ContourTreeAugmented& filter, bool dataFieldIsSorted, const vtkm::cont::UnknownArrayHandle& arr);
30+
void distributed_analysis(const vtkm::cont::PartitionedDataSet& result, int rank);
3031

3132
std::string m_field_name;
3233
int m_levels;

‎src/tests/vtkh/t_vtk-h_contour_tree_par.cpp‎

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ DEFINE_MPI_CAST(MPI_Comm)
4343

4444
#endif
4545

46-
using ValueType = vtkm::Float64;
46+
using ValueType = vtkm::Float32;
4747

4848
inline vtkm::IdComponent FindSplitAxis(vtkm::Id3 globalSize)
4949
{
@@ -446,11 +446,19 @@ TEST(vtkh_contour_tree, vtkh_contour_tree)
446446
std::vector<double> isoValues = marcher.GetIsoValues();
447447
std::sort(isoValues.begin(), isoValues.end());
448448

449-
EXPECT_FLOAT_EQ(isoValues[0], 1e-05);
450-
EXPECT_FLOAT_EQ(isoValues[1], 82);
451-
EXPECT_FLOAT_EQ(isoValues[2], 133);
452-
EXPECT_FLOAT_EQ(isoValues[3], 168);
453-
EXPECT_FLOAT_EQ(isoValues[4], 177);
449+
if( mpiSize == 1 ) {
450+
EXPECT_FLOAT_EQ(isoValues[0], 1e-05);
451+
EXPECT_FLOAT_EQ(isoValues[1], 82);
452+
EXPECT_FLOAT_EQ(isoValues[2], 133);
453+
EXPECT_FLOAT_EQ(isoValues[3], 168);
454+
EXPECT_FLOAT_EQ(isoValues[4], 177);
455+
} else {
456+
EXPECT_FLOAT_EQ(isoValues[0], 1e-05);
457+
EXPECT_FLOAT_EQ(isoValues[1], 17.00001);
458+
EXPECT_FLOAT_EQ(isoValues[2], 82);
459+
EXPECT_FLOAT_EQ(isoValues[3], 133);
460+
EXPECT_FLOAT_EQ(isoValues[4], 176);
461+
}
454462

455463
vtkh::DataSet *output = marcher.GetOutput();
456464
if( output )

0 commit comments

Comments
 (0)