Refactor merging to be restricted to only common vertices

This commit is contained in:
Gunther H. Weber 2021-06-16 16:00:29 -07:00
parent 47cc80e42c
commit 32dcd3d0c5
4 changed files with 369 additions and 181 deletions

@ -70,6 +70,7 @@
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ArrayPortalToIterators.h>
#include <vtkm/cont/ArrayRangeCompute.h>
@ -90,9 +91,10 @@
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedArrayWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedNeighborsWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyNeighborsToPackedArray.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindDuplicateInOtherWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/GetArcFromDecorator.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeIntoCombinedNeighborsWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeSortedListsWithoutDuplicatesWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
@ -414,6 +416,35 @@ inline ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
#endif
}
template <typename PT1, typename PT2, typename PT3, typename PT4>
inline void CopyArrayByIndices(const PT1& srcArray,
const PT2& srcIndices,
PT3& tgtArray,
const PT4& tgtIndices)
{
VTKM_ASSERT(srcIndices.GetNumberOfValues() == tgtIndices.GetNumberOfValues());
auto srcPermutation = make_ArrayHandlePermutation(srcIndices, srcArray);
auto tgtPermuation = make_ArrayHandlePermutation(tgtIndices, tgtArray);
vtkm::cont::Algorithm::Copy(srcPermutation, tgtPermuation);
}
template <typename PT1, typename PT2, typename PT3, typename PT4>
inline void CopyVecArrayByIndices(const PT1& srcArray,
const PT2& srcIndices,
PT3& tgtArray,
const PT4& tgtIndices)
{
VTKM_ASSERT(srcIndices.GetNumberOfValues() == tgtIndices.GetNumberOfValues());
auto srcPermutation = make_ArrayHandlePermutation(srcIndices, srcArray);
auto tgtPermuation = make_ArrayHandlePermutation(tgtIndices, tgtArray);
vtkm::cont::Invoker invoke;
invoke(
contourtree_mesh_inc_ns::CopyIntoCombinedNeighborsWorklet{}, srcPermutation, tgtPermuation);
// Why doesn't the following copy work instead?
//vtkm::cont::Algorithm::Copy(srcPermutation, tgtPermuation);
}
// Initalize the contour tree from the arcs array and sort order
template <typename FieldType>
@ -657,125 +688,195 @@ inline void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& ot
<< ": " << timer.GetElapsedTime() << " seconds" << std::endl;
timer.Start();
// Combined neighbor counts (including potential duplicates)
IdArrayType combinedNeighborCounts;
vtkm::cont::Algorithm::Fill(combinedNeighborCounts, vtkm::Id{ 0 }, numVerticesCombined);
// First, copy our neighbor counts into the correct locations
auto neighborCountsThis = make_ArrayHandleOffsetsToNumComponents(this->NeighborOffsets);
auto permutedCombinedNeighborCountsThis =
vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, combinedNeighborCounts);
vtkm::cont::Algorithm::Copy(neighborCountsThis, permutedCombinedNeighborCountsThis);
// Now, add the neighbor counts from the other mesh at the correct positions
auto neighborCountsOther = make_ArrayHandleOffsetsToNumComponents(other.NeighborOffsets);
auto permutedCombinedNeighborCountsOther =
vtkm::cont::make_ArrayHandlePermutation(otherToCombinedSortOrder, combinedNeighborCounts);
this->Invoke(contourtree_mesh_inc_ns::AddToArrayElementsWorklet{},
permutedCombinedNeighborCountsOther,
neighborCountsOther);
timingsStream << " " << std::setw(38) << std::left << "Create CombinedOtherStartIndex"
<< ": " << timer.GetElapsedTime() << " seconds" << std::endl;
timer.Start();
// Compute offsets from counts
IdArrayType combinedNeighborOffsets;
vtkm::Id nCombinedNeighbours;
vtkm::cont::ConvertNumComponentsToOffsets(
combinedNeighborCounts, combinedNeighborOffsets, nCombinedNeighbours);
//PrintIndices("combinedNeighborCounts", combinedNeighborCounts);
//PrintIndices("combinedNeighborOffsets", combinedNeighborOffsets);
// Create combined neighbor connectivity array and grouped view on it
IdArrayType combinedNeighborConnectivity;
combinedNeighborConnectivity.Allocate(nCombinedNeighbours);
auto combinedNeighborGroups = vtkm::cont::make_ArrayHandleGroupVecVariable(
combinedNeighborConnectivity, combinedNeighborOffsets);
//printNumComponents(combinedNeighborGroups);
// Create array to hold merged neighbor counts (after eliminating duplicate neigbors)
vtkm::cont::ArrayHandle<vtkm::IdComponent> mergedNeighborCounts;
vtkm::cont::Algorithm::Fill(mergedNeighborCounts, vtkm::IdComponent{ 0 }, numVerticesCombined);
// Copy connectivity array from this mesh into combined neighbor connectivity
// The following decorator converts vertiex indices from our mesh to the combined mesh
auto thisNeighborConnectivityGlobal =
// Map neighbor IDs to global ID (ID in the combined) and group them
auto neighborConnectivityGlobalThis =
make_ArrayHandleDecorator(this->NeighborConnectivity.GetNumberOfValues(),
mesh_dem_contourtree_mesh_inc::ApplyLookupTableDecoratorImpl{},
this->NeighborConnectivity,
thisToCombinedSortOrder);
auto neighborConnectivityGlobalGroupsThis = vtkm::cont::make_ArrayHandleGroupVecVariable(
neighborConnectivityGlobalThis, this->NeighborOffsets);
//PrintArray("thisNeighborConnectivityGlobal", thisNeighborConnectivityGlobal);
//PrintArray("this->NeighborOffsets", this->NeighborOffsets);
auto thisNeighborConnectivityGlobalGroups = vtkm::cont::make_ArrayHandleGroupVecVariable(
thisNeighborConnectivityGlobal, this->NeighborOffsets);
//printNumComponents(thisNeighborConnectivityGlobalGroups);
auto permutedMergedNeighborCountsThis =
vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, mergedNeighborCounts);
auto permutedCombinedNeighborGroupsThis =
vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, combinedNeighborGroups);
this->Invoke(contourtree_mesh_inc_ns::CopyIntoCombinedNeighborsWorklet{},
thisNeighborConnectivityGlobalGroups,
permutedCombinedNeighborGroupsThis,
permutedMergedNeighborCountsThis);
//std::cout << "After CopyIntoCombinedNeighborsWorklet" << std::endl;
//PrintArray("combinedNeighborConnectivity", combinedNeighborConnectivity);
//PrintArray("combinedNeighborOffsets", combinedNeighborOffsets);
//printNumComponents(combinedNeighborGroups);
// Insert neighbors from other mesh into combined neighbor connectivity (without duplicates)
// The following decorator converts vertiex indices from our mesh to the combined mesh
//
auto otherNeighborConnectivityGlobal =
auto neighborConnectivityGlobalOther =
make_ArrayHandleDecorator(other.NeighborConnectivity.GetNumberOfValues(),
mesh_dem_contourtree_mesh_inc::ApplyLookupTableDecoratorImpl{},
other.NeighborConnectivity,
otherToCombinedSortOrder);
auto neighborConnectivityGlobalGroupsOther = vtkm::cont::make_ArrayHandleGroupVecVariable(
neighborConnectivityGlobalOther, other.NeighborOffsets);
//PrintArray("otherNeighborConnectivityGlobal", otherNeighborConnectivityGlobal);
//PrintArray("other.NeighborOffsets", other.NeighborOffsets);
auto otherNeighborConnectivityGlobalGroups = vtkm::cont::make_ArrayHandleGroupVecVariable(
otherNeighborConnectivityGlobal, other.NeighborOffsets);
//printNumComponents(otherNeighborConnectivityGlobalGroups);
auto permutedMergedNeighborCountsOther =
vtkm::cont::make_ArrayHandlePermutation(otherToCombinedSortOrder, mergedNeighborCounts);
auto permutedCombinedNeighborGroupsOther =
vtkm::cont::make_ArrayHandlePermutation(otherToCombinedSortOrder, combinedNeighborGroups);
this->Invoke(contourtree_mesh_inc_ns::MergeIntoCombinedNeighborsWorklet{},
otherNeighborConnectivityGlobalGroups,
permutedCombinedNeighborGroupsOther,
permutedMergedNeighborCountsOther);
// Merge the two neighborhood connecitivy lists. First, we split neighbor connectivity
// into three groups (i) vertices only in this, (ii) vertices only in other, (iii)
// vertices in both meshes. We then compute cobmined neighbor connectivity for vertices
// in both meshes. Finally, we copy them into the combined array.
//std::cout << "After MergeIntoCombinedNeighborsWorklet" << std::endl;
//PrintArray("combinedNeighborConnectivity", combinedNeighborConnectivity);
//PrintArray("combinedNeighborOffsets", combinedNeighborOffsets);
//printNumComponents(combinedNeighborGroups);
timingsStream << " " << std::setw(38) << std::left << "Update CombinedNeighbours"
// Split vertices into groups (i) uniuqe this, (ii) unique other, (iii) in both
// ... compute arrays containing 1 if the element is in the other respective array
vtkm::cont::ArrayHandle<vtkm::IdComponent> thisToCombinedSortOrderIsDuplicate;
thisToCombinedSortOrderIsDuplicate.Allocate(thisToCombinedSortOrder.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::IdComponent> otherToCombinedSortOrderIsDuplicate;
vtkm::cont::Algorithm::Fill(otherToCombinedSortOrderIsDuplicate,
vtkm::IdComponent{ 0 },
otherToCombinedSortOrder.GetNumberOfValues());
this->Invoke(contourtree_mesh_inc_ns::FindDuplicateInOtherWorklet{},
thisToCombinedSortOrder,
otherToCombinedSortOrder,
thisToCombinedSortOrderIsDuplicate,
otherToCombinedSortOrderIsDuplicate);
#ifdef DEBUG_PRINT
PrintIndices("thisToCombinedSortOrderIsDuplicate", thisToCombinedSortOrderIsDuplicate);
PrintIndices("otherToCombinedSortOrderIsDuplicate", otherToCombinedSortOrderIsDuplicate);
#endif
// ... create lists for all groups to be used to restrict operations to them
vtkm::cont::ArrayHandleIndex indicesThis(thisToCombinedSortOrder.GetNumberOfValues());
vtkm::cont::ArrayHandleIndex indicesOther(otherToCombinedSortOrder.GetNumberOfValues());
// ... helper function, basically negates criterion for CopyIf
struct IsUnique
{
VTKM_EXEC_CONT bool operator()(vtkm::IdComponent isInOther) { return isInOther == 0; }
};
IdArrayType indicesThisUnique, indicesThisDuplicate;
vtkm::cont::Algorithm::CopyIf(
indicesThis, thisToCombinedSortOrderIsDuplicate, indicesThisUnique, IsUnique{});
vtkm::cont::Algorithm::CopyIf(
indicesThis, thisToCombinedSortOrderIsDuplicate, indicesThisDuplicate);
#ifdef DEBUG_PRINT
PrintIndices("indicesThisUnique", indicesThisUnique);
PrintIndices("indicesThisDuplicate", indicesThisDuplicate);
#endif
IdArrayType indicesOtherUnique, indicesOtherDuplicate;
vtkm::cont::Algorithm::CopyIf(
indicesOther, otherToCombinedSortOrderIsDuplicate, indicesOtherUnique, IsUnique{});
vtkm::cont::Algorithm::CopyIf(
indicesOther, otherToCombinedSortOrderIsDuplicate, indicesOtherDuplicate);
#ifdef DEBUG_PRINT
PrintIndices("indicesOtherUnique", indicesOtherUnique);
PrintIndices("indicesOtherDuplicate", indicesOtherDuplicate);
#endif
VTKM_ASSERT(indicesThisDuplicate.GetNumberOfValues() ==
indicesOtherDuplicate.GetNumberOfValues());
// Merge the neighbor groups for vertices that occur in both meshes
// ... compute combined counts (with duplicates)
auto neighborCountsThis = make_ArrayHandleOffsetsToNumComponents(this->NeighborOffsets);
auto permutedNeighborCountsThis =
vtkm::cont::make_ArrayHandlePermutation(indicesThisDuplicate, neighborCountsThis);
auto neighborCountsOther = make_ArrayHandleOffsetsToNumComponents(other.NeighborOffsets);
auto permutedNeighborCountsOther =
vtkm::cont::make_ArrayHandlePermutation(indicesOtherDuplicate, neighborCountsOther);
vtkm::cont::ArrayHandle<vtkm::IdComponent> combinedCommonNeighborCountSums;
vtkm::cont::Algorithm::Transform(permutedNeighborCountsThis,
permutedNeighborCountsOther,
combinedCommonNeighborCountSums,
vtkm::Sum());
// ... merge sorted lists
// ...... create output arrays/groups
vtkm::Id unpackedCombinedCommonNeighborConnectivitySize;
IdArrayType unpackedCombinedCommonNeighborOffsets;
vtkm::cont::ConvertNumComponentsToOffsets(combinedCommonNeighborCountSums,
unpackedCombinedCommonNeighborOffsets,
unpackedCombinedCommonNeighborConnectivitySize);
IdArrayType unpackedCombinedCommonNeighborConnectivity;
unpackedCombinedCommonNeighborConnectivity.Allocate(
unpackedCombinedCommonNeighborConnectivitySize);
auto unpackedCombinedCommonNeighborConnectivityGroups = make_ArrayHandleGroupVecVariable(
unpackedCombinedCommonNeighborConnectivity, unpackedCombinedCommonNeighborOffsets);
// ....... create permuted input arrays/groups
auto permutedNeighborConnectivityGlobalGroupsThis =
make_ArrayHandlePermutation(indicesThisDuplicate, neighborConnectivityGlobalGroupsThis);
auto permutedNeighborConnectivityGlobalGroupsOther =
make_ArrayHandlePermutation(indicesOtherDuplicate, neighborConnectivityGlobalGroupsOther);
// ........ create array for actual counts of unique neighbors
vtkm::cont::ArrayHandle<vtkm::IdComponent> packedCombinedCommonNeighborCounts;
packedCombinedCommonNeighborCounts.Allocate(combinedCommonNeighborCountSums.GetNumberOfValues());
// ........ perform merge
this->Invoke(contourtree_mesh_inc_ns::MergeSortedListsWithoutDuplicatesWorklet{},
permutedNeighborConnectivityGlobalGroupsThis,
permutedNeighborConnectivityGlobalGroupsOther,
unpackedCombinedCommonNeighborConnectivityGroups,
packedCombinedCommonNeighborCounts);
// ... pack sorted lists
// ...... create the new offsets array for the merged groups (without duplicates).
vtkm::Id packedCombinedCommonNeighborConnectivitySize;
vtkm::cont::ArrayHandle<vtkm::Id> packedCombinedCommonNeighborOffsets;
vtkm::cont::ConvertNumComponentsToOffsets(packedCombinedCommonNeighborCounts,
packedCombinedCommonNeighborOffsets,
packedCombinedCommonNeighborConnectivitySize);
// ...... create a new grouped array for the packed connectivity
vtkm::cont::ArrayHandle<vtkm::Id> packedCombinedCommonNeighborConnectivity;
packedCombinedCommonNeighborConnectivity.Allocate(packedCombinedCommonNeighborConnectivitySize);
auto packedCommonNeighborConnectivityGroups = vtkm::cont::make_ArrayHandleGroupVecVariable(
packedCombinedCommonNeighborConnectivity, packedCombinedCommonNeighborOffsets);
// ...... copy data to the packed array.
this->Invoke(contourtree_mesh_inc_ns::CopyNeighborsToPackedArray{},
unpackedCombinedCommonNeighborConnectivityGroups,
packedCommonNeighborConnectivityGroups);
// Create array for all three groups
// ... create combined counts array
IdArrayType combinedNeighborCounts;
combinedNeighborCounts.Allocate(numVerticesCombined);
auto thisOnlyToCombinedSortOrder =
make_ArrayHandlePermutation(indicesThisUnique, thisToCombinedSortOrder);
auto otherOnlyToCombinedSortOrder =
make_ArrayHandlePermutation(indicesOtherUnique, otherToCombinedSortOrder);
auto commonToCombinedSortOrder =
make_ArrayHandlePermutation(indicesThisDuplicate, thisToCombinedSortOrder);
CopyArrayByIndices(
neighborCountsThis, indicesThisUnique, combinedNeighborCounts, thisOnlyToCombinedSortOrder);
CopyArrayByIndices(
neighborCountsOther, indicesOtherUnique, combinedNeighborCounts, otherOnlyToCombinedSortOrder);
auto commonCombinedNeighborCounts =
make_ArrayHandlePermutation(commonToCombinedSortOrder, combinedNeighborCounts);
vtkm::cont::Algorithm::Copy(packedCombinedCommonNeighborCounts, commonCombinedNeighborCounts);
// ... create offsets and allocate combinedNeighborConnectivity array
vtkm::Id combinedNeighborConnectivitySize;
vtkm::cont::ArrayHandle<vtkm::Id> combinedNeighborOffsets;
vtkm::cont::ConvertNumComponentsToOffsets(
combinedNeighborCounts, combinedNeighborOffsets, combinedNeighborConnectivitySize);
IdArrayType combinedNeighborConnectivity;
combinedNeighborConnectivity.Allocate(combinedNeighborConnectivitySize);
auto combinedNeighborConnectivityGroups =
make_ArrayHandleGroupVecVariable(combinedNeighborConnectivity, combinedNeighborOffsets);
// ... copy the connectivity data including previously merged lists
CopyVecArrayByIndices(neighborConnectivityGlobalGroupsThis,
indicesThisUnique,
combinedNeighborConnectivityGroups,
thisOnlyToCombinedSortOrder);
CopyVecArrayByIndices(neighborConnectivityGlobalGroupsOther,
indicesOtherUnique,
combinedNeighborConnectivityGroups,
otherOnlyToCombinedSortOrder);
auto commonCombinedNeighborConnectivityGroups =
make_ArrayHandlePermutation(commonToCombinedSortOrder, combinedNeighborConnectivityGroups);
this->Invoke(contourtree_mesh_inc_ns::CopyIntoCombinedNeighborsWorklet{},
packedCommonNeighborConnectivityGroups,
commonCombinedNeighborConnectivityGroups);
// Why doesn't the following copy work instead?
// vtkm::cont::Algorithm::Copy(packedCommonNeighborConnectivityGroups, commonCombinedNeighborConnectivityGroups);
timingsStream << " " << std::setw(38) << std::left << "Compute CombinedNeighborConnectivity"
<< ": " << timer.GetElapsedTime() << " seconds" << std::endl;
timer.Start();
// Create the new offsets array for the merged groups (without duplicates).
vtkm::Id packedCombinedNeighborConnectivitySize;
vtkm::cont::ArrayHandle<vtkm::Id> packedCombinedNeighborOffsets;
vtkm::cont::ConvertNumComponentsToOffsets(
mergedNeighborCounts, packedCombinedNeighborOffsets, packedCombinedNeighborConnectivitySize);
// Create a new grouped array for the packed connectivity
// Note that you have to pre-allocate the connectivity array.
vtkm::cont::ArrayHandle<vtkm::Id> packedCombinedNeighborConnectivity;
packedCombinedNeighborConnectivity.Allocate(packedCombinedNeighborConnectivitySize);
auto packedNeighborGroups = vtkm::cont::make_ArrayHandleGroupVecVariable(
packedCombinedNeighborConnectivity, packedCombinedNeighborOffsets);
// Copy data to the packed array.
this->Invoke(contourtree_mesh_inc_ns::CopyNeighborsToPackedArray{},
combinedNeighborGroups,
packedNeighborGroups);
// Compute combined global mesh index arrays
IdArrayType combinedGlobalMeshIndex;
combinedGlobalMeshIndex.Allocate(numVerticesCombined);
@ -815,8 +916,8 @@ inline void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& ot
// Swap in combined version. VTKM ArrayHandles are smart so we can just swap in the new for the old
this->SortedValues = combinedSortedValues;
this->GlobalMeshIndex = combinedGlobalMeshIndex;
this->NeighborConnectivity = packedCombinedNeighborConnectivity;
this->NeighborOffsets = packedCombinedNeighborOffsets;
this->NeighborConnectivity = combinedNeighborConnectivity;
this->NeighborOffsets = combinedNeighborOffsets;
this->NumVertices = SortedValues.GetNumberOfValues();
this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);

@ -77,21 +77,18 @@ namespace mesh_dem_contourtree_mesh_inc
struct CopyIntoCombinedNeighborsWorklet : public vtkm::worklet::WorkletMapField
{
// TODO/FIXME: Even though the second parameter is only used for output,
// GetNumberOfComponents() in the assert only returns the correct value
// if defined as FieldInOut. Possibly skip assert, but then the later
// assign in the for loop may fail. Check with Ken Moreland what the
// correct way to implement this is.
using ControlSignature = void(FieldIn, FieldInOut, FieldOut);
// NOTE: Even though the second parameter is only used for output, it needs
// to be declared as FieldInOut as we use an ArrayHandleGroupVecVariable for
// output, which breaks some assumptions. Otherwise GetNumberOfComponents()
// and indexing do not work properly (see issue #626) resulting in a crash.
using ControlSignature = void(FieldIn, FieldInOut);
template <typename InGroupType, typename OutGroupType>
VTKM_EXEC void operator()(const InGroupType& newNeighbors,
OutGroupType& combinedNeighbors,
vtkm::IdComponent& actualGroupSize) const
VTKM_EXEC void operator()(const InGroupType& newNeighbors, OutGroupType& combinedNeighbors) const
{
actualGroupSize = newNeighbors.GetNumberOfComponents();
VTKM_ASSERT(actualGroupSize <= combinedNeighbors.GetNumberOfComponents());
for (vtkm::IdComponent index = 0; index < actualGroupSize; ++index)
vtkm::IdComponent groupSize = newNeighbors.GetNumberOfComponents();
VTKM_ASSERT(groupSize == combinedNeighbors.GetNumberOfComponents());
for (vtkm::IdComponent index = 0; index < groupSize; ++index)
{
combinedNeighbors[index] = newNeighbors[index];
}

@ -0,0 +1,114 @@
//============================================================================
// 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 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
// 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)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_find_duplicate_in_other_worklet_h
#define vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_find_duplicate_in_other_worklet_h
#include <vtkm/LowerBound.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace mesh_dem_contourtree_mesh_inc
{
struct FindDuplicateInOtherWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn, WholeArrayIn, FieldOut, WholeArrayOut);
template <typename InputType,
typename InputArrayPortalType,
typename OutputType,
typename OutputArrayPortalType>
VTKM_EXEC void operator()(const InputType& value,
const InputArrayPortalType& otherArrayPortal,
OutputType& isDuplicate,
OutputArrayPortalType& otherDuplicatePortal) const
{
vtkm::Id posInOther = vtkm::LowerBound(otherArrayPortal, value);
if (posInOther < otherArrayPortal.GetNumberOfValues() &&
otherArrayPortal.Get(posInOther) == value)
{
isDuplicate = 1;
otherDuplicatePortal.Set(posInOther, 1);
}
else
{
isDuplicate = 0;
}
}
}; // FindDuplicateInOtherWorklet
} // namespace mesh_dem_contourtree_mesh_inc
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm
#endif

@ -60,8 +60,8 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_add_to_array_elements__worklet_worklet_h
#define vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_add_to_array_elements__worklet_worklet_h
#ifndef vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_merge_sorted_lists_worklet_h
#define vtk_m_worklet_contourtree_augmented_contourtree_mesh_inc_merge_sorted_lists_worklet_h
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
@ -75,75 +75,51 @@ namespace contourtree_augmented
namespace mesh_dem_contourtree_mesh_inc
{
struct MergeIntoCombinedNeighborsWorklet : public vtkm::worklet::WorkletMapField
struct MergeSortedListsWithoutDuplicatesWorklet : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn, FieldInOut, FieldInOut);
using ControlSignature = void(FieldIn, FieldIn, FieldOut, FieldOut);
template <typename InGroupType, typename InOutGroupType>
VTKM_EXEC void operator()(const InGroupType& newNeighbors,
InOutGroupType& combinedNeighbors,
vtkm::IdComponent& actualGroupSize) const
template <typename InGroupType, typename OutGroupType>
VTKM_EXEC void operator()(const InGroupType& list1,
const InGroupType& list2,
OutGroupType& combinedList,
vtkm::IdComponent& numberOfUniqueElements) const
{
VTKM_ASSERT(actualGroupSize + newNeighbors.GetNumberOfComponents() <=
combinedNeighbors.GetNumberOfComponents());
VTKM_ASSERT(list1.GetNumberOfComponents() + list2.GetNumberOfComponents() <=
combinedList.GetNumberOfComponents());
if (actualGroupSize == 0)
numberOfUniqueElements = 0;
vtkm::IdComponent pos1 = 0;
vtkm::IdComponent pos2 = 0;
while (pos1 < list1.GetNumberOfComponents() && pos2 < list2.GetNumberOfComponents())
{
// Most common case, no neighbor list present, yet -> just copy
actualGroupSize = newNeighbors.GetNumberOfComponents();
for (vtkm::IdComponent index = 0; index < actualGroupSize; ++index)
if (list1[pos1] < list2[pos2])
{
combinedNeighbors[index] = newNeighbors[index];
combinedList[numberOfUniqueElements++] = list1[pos1++];
}
else if (list1[pos1] == list2[pos2])
{
combinedList[numberOfUniqueElements++] = list1[pos1++];
++pos2;
}
else
{
// Shared vertex -> Need to merge sorted list of new neighbors into sorted
// list of already existing neighbors
// TODO/FIXME: Better way to merge two sorted lists?
vtkm::IdComponent numToInsert = newNeighbors.GetNumberOfComponents();
// Define insert pos here. Since both lists are sorted, subsequent inserts
// will always occur after the previous one. So it makes sense to start
// the search for a new insert position at the last one.
vtkm::IdComponent insertPos = 0;
for (vtkm::IdComponent idxToInsert = 0; idxToInsert < numToInsert; ++idxToInsert)
VTKM_ASSERT(list1[pos1] > list2[pos2]);
combinedList[numberOfUniqueElements++] = list2[pos2++];
}
}
// Either list1 or list2 may have additional elements but not both, so the following is safe
while (pos1 < list1.GetNumberOfComponents())
{
//std::cout << "actualGroupSize = " << actualGroupSize;
//std::cout << " inserting " << newNeighbors[idxToInsert] << std::endl;
//std::cout << insertPos << " " << (insertPos < actualGroupSize) << " "
// << (combinedNeighbors[insertPos] < newNeighbors[idxToInsert]) << std::endl;
while (insertPos < actualGroupSize &&
combinedNeighbors[insertPos] < newNeighbors[idxToInsert])
combinedList[numberOfUniqueElements++] = list1[pos1++];
}
while (pos2 < list2.GetNumberOfComponents())
{
//std::cout << insertPos << " " << (insertPos < actualGroupSize) << " "
// << (combinedNeighbors[insertPos] < newNeighbors[idxToInsert]) << std::endl;
++insertPos;
}
//std::cout << "Insert pos is " << insertPos << std::endl;
if ((insertPos >= actualGroupSize) ||
combinedNeighbors[insertPos] != newNeighbors[idxToInsert])
{
// Only insert non-duplicate elements
// Shift elements one back
for (vtkm::IdComponent idx = actualGroupSize - 1; idx >= insertPos; --idx)
{
//std::cout << "Shifting " << combinedNeighbors[idx] << " from " << idx << " to "
// << idx + 1 << std::endl;
combinedNeighbors[idx + 1] = combinedNeighbors[idx];
}
//std::cout << "Saving " << newNeighbors[idxToInsert] << " to pos " << insertPos
// << std::endl;
// Insert element and update group size
combinedNeighbors[insertPos] = newNeighbors[idxToInsert];
actualGroupSize += 1;
// Since both lists are sorted, the lowest possible insert position is behind the current element
insertPos += 1;
combinedList[numberOfUniqueElements++] = list2[pos2++];
}
}
}
}
}; // MergeIntoCombinedNeighborsWorklet
}; // MergeSortedListsWithoutDuplicatesWorklet
} // namespace mesh_dem_contourtree_mesh_inc
} // namespace contourtree_augmented