diff --git a/vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h b/vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h index 9168bb982..7c43d5cb6 100644 --- a/vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h +++ b/vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h @@ -70,6 +70,7 @@ #include #include #include +#include #include #include #include @@ -90,9 +91,10 @@ #include #include #include +#include #include #include -#include +#include #include #include #include @@ -114,7 +116,7 @@ VTKM_THIRDPARTY_POST_INCLUDE namespace contourtree_mesh_inc_ns = vtkm::worklet::contourtree_augmented::mesh_dem_contourtree_mesh_inc; -//#define DEBUG_PRINT +// #define DEBUG_PRINT namespace vtkm { @@ -414,6 +416,35 @@ inline ContourTreeMesh::ContourTreeMesh(const IdArrayType& nodes, #endif } +template +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 +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 @@ -657,125 +688,195 @@ inline void ContourTreeMesh::MergeWith(ContourTreeMesh& 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 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 thisToCombinedSortOrderIsDuplicate; + thisToCombinedSortOrderIsDuplicate.Allocate(thisToCombinedSortOrder.GetNumberOfValues()); + vtkm::cont::ArrayHandle 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 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 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 packedCombinedCommonNeighborOffsets; + vtkm::cont::ConvertNumComponentsToOffsets(packedCombinedCommonNeighborCounts, + packedCombinedCommonNeighborOffsets, + packedCombinedCommonNeighborConnectivitySize); + + // ...... create a new grouped array for the packed connectivity + vtkm::cont::ArrayHandle 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 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 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 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::MergeWith(ContourTreeMesh& 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); diff --git a/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedNeighborsWorklet.h b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedNeighborsWorklet.h index df2e38a99..dda125a3b 100644 --- a/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedNeighborsWorklet.h +++ b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CopyIntoCombinedNeighborsWorklet.h @@ -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 - 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]; } diff --git a/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindDuplicateInOtherWorklet.h b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindDuplicateInOtherWorklet.h new file mode 100644 index 000000000..573335516 --- /dev/null +++ b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindDuplicateInOtherWorklet.h @@ -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 +#include +#include + +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 + + 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 diff --git a/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeIntoCombinedNeighborsWorklet.h b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeSortedListsWithoutDuplicatesWorklet.h similarity index 54% rename from vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeIntoCombinedNeighborsWorklet.h rename to vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeSortedListsWithoutDuplicatesWorklet.h index 3190c7b4b..77115928e 100644 --- a/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeIntoCombinedNeighborsWorklet.h +++ b/vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeSortedListsWithoutDuplicatesWorklet.h @@ -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 #include @@ -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 - VTKM_EXEC void operator()(const InGroupType& newNeighbors, - InOutGroupType& combinedNeighbors, - vtkm::IdComponent& actualGroupSize) const + template + 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 + { + VTKM_ASSERT(list1[pos1] > list2[pos2]); + combinedList[numberOfUniqueElements++] = list2[pos2++]; } } - else + // Either list1 or list2 may have additional elements but not both, so the following is safe + while (pos1 < list1.GetNumberOfComponents()) { - // 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) - { - //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]) - { - //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++] = list1[pos1++]; + } + while (pos2 < list2.GetNumberOfComponents()) + { + combinedList[numberOfUniqueElements++] = list2[pos2++]; } } -}; // MergeIntoCombinedNeighborsWorklet - +}; // MergeSortedListsWithoutDuplicatesWorklet } // namespace mesh_dem_contourtree_mesh_inc } // namespace contourtree_augmented