Commit 2714da8f authored by Abhishek Yenpure's avatar Abhishek Yenpure
Browse files

Adding changes suggested by Rob on GitLab

- Seperate BoundingIntervalHierarchy.h into
  - BoundingIntervalHierarchy.h
  - BoundingIntervalHierarchy.hxx
  - BoundingIntervalHierarchy.cxx
parent ab81c633
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/BoundingIntervalHierarchy.h>
#include <vtkm/cont/BoundingIntervalHierarchy.hxx>
namespace vtkm
{
namespace cont
{
VTKM_CONT
void BoundingIntervalHierarchy::Build()
{
BuildFunctor functor(this);
vtkm::cont::TryExecute(functor);
}
VTKM_CONT
const vtkm::exec::CellLocator* BoundingIntervalHierarchy::PrepareForExecutionImpl(
const vtkm::Int8 device) const
{
using DeviceList = VTKM_DEFAULT_DEVICE_ADAPTER_LIST_TAG;
const vtkm::exec::CellLocator* toReturn;
vtkm::cont::internal::FindDeviceAdapterTagAndCall(
device, DeviceList(), PrepareForExecutionFunctor(), this, &toReturn);
return toReturn;
}
}
}
This diff is collapsed.
This diff is collapsed.
......@@ -103,6 +103,7 @@ set(headers
set(template_sources
ArrayHandle.hxx
ArrayRangeCompute.hxx
BoundingIntervalHierarchy.hxx
CellSetExplicit.hxx
CellSetStructured.hxx
CoordinateSystem.hxx
......@@ -114,6 +115,7 @@ set(template_sources
set(sources
ArrayHandle.cxx
AssignerMultiBlock.cxx
BoundingIntervalHierarchy.cxx
BoundsCompute.cxx
BoundsGlobalCompute.cxx
CellSet.cxx
......
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_worklet_spatialstructure_BoundingIntervalHierarchy_h
#define vtk_m_worklet_spatialstructure_BoundingIntervalHierarchy_h
#include <type_traits>
#include <vtkm/Bounds.h>
#include <vtkm/Types.h>
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ArrayHandleReverse.h>
#include <vtkm/cont/ArrayHandleTransform.h>
#include <vtkm/cont/BoundingIntervalHierarchyNode.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/cuda/DeviceAdapterCuda.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
namespace vtkm
{
namespace worklet
{
namespace spatialstructure
{
struct TreeNode
{
vtkm::FloatDefault LMax;
vtkm::FloatDefault RMin;
vtkm::IdComponent Dimension;
TreeNode()
: LMax()
, RMin()
, Dimension()
{
}
}; // struct TreeNode
struct SplitProperties
{
vtkm::FloatDefault Plane;
vtkm::Id NumLeftPoints;
vtkm::Id NumRightPoints;
vtkm::FloatDefault LMax;
vtkm::FloatDefault RMin;
vtkm::FloatDefault Cost;
SplitProperties()
: Plane()
, NumLeftPoints()
, NumRightPoints()
, LMax()
, RMin()
, Cost()
{
}
}; // struct SplitProperties
struct CellRangesExtracter : public vtkm::worklet::WorkletMapPointToCell
{
typedef void ControlSignature(CellSetIn,
WholeArrayIn<>,
FieldOutCell<>,
FieldOutCell<>,
FieldOutCell<>,
FieldOutCell<>,
FieldOutCell<>,
FieldOutCell<>);
typedef void ExecutionSignature(_1, PointIndices, _2, _3, _4, _5, _6, _7, _8);
template <typename CellShape, typename PointIndicesVec, typename PointsPortal>
VTKM_EXEC void operator()(CellShape vtkmNotUsed(shape),
const PointIndicesVec& pointIndices,
const PointsPortal& points,
vtkm::Range& rangeX,
vtkm::Range& rangeY,
vtkm::Range& rangeZ,
vtkm::FloatDefault& centerX,
vtkm::FloatDefault& centerY,
vtkm::FloatDefault& centerZ) const
{
vtkm::Bounds bounds;
vtkm::VecFromPortalPermute<PointIndicesVec, PointsPortal> cellPoints(&pointIndices, points);
vtkm::IdComponent numPoints = cellPoints.GetNumberOfComponents();
for (vtkm::IdComponent i = 0; i < numPoints; ++i)
{
bounds.Include(cellPoints[i]);
}
rangeX = bounds.X;
rangeY = bounds.Y;
rangeZ = bounds.Z;
vtkm::Vec<vtkm::FloatDefault, 3> center = bounds.Center();
centerX = center[0];
centerY = center[1];
centerZ = center[2];
}
}; // struct CellRangesExtracter
struct LEQWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldOut<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::FloatDefault& value,
const vtkm::FloatDefault& planeValue,
vtkm::Id& leq,
vtkm::Id& r) const
{
leq = value <= planeValue;
r = !leq;
}
}; // struct LEQWorklet
template <bool LEQ>
struct FilterRanges : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::FloatDefault& value,
const vtkm::FloatDefault& planeValue,
const vtkm::Range& cellBounds,
vtkm::Range& outBounds) const
{
if (LEQ)
{
outBounds = (value <= planeValue) ? cellBounds : vtkm::Range();
}
else
{
outBounds = (value > planeValue) ? cellBounds : vtkm::Range();
}
}
}; // struct FilterRanges
struct SplitPlaneCalculatorWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2);
using InputDomain = _1;
VTKM_CONT
SplitPlaneCalculatorWorklet(vtkm::IdComponent planeIdx, vtkm::IdComponent numPlanes)
: Scale(static_cast<vtkm::FloatDefault>(planeIdx + 1) /
static_cast<vtkm::FloatDefault>(numPlanes + 1))
{
}
VTKM_EXEC
void operator()(const vtkm::Range& range, vtkm::FloatDefault& splitPlane) const
{
splitPlane = range.Min + Scale * (range.Max - range.Min);
}
vtkm::FloatDefault Scale;
};
struct SplitPropertiesCalculator : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldIn<>,
FieldIn<>,
FieldIn<>,
WholeArrayInOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, InputIndex);
using InputDomain = _1;
VTKM_CONT
SplitPropertiesCalculator(vtkm::IdComponent index, vtkm::Id stride)
: Index(index)
, Stride(stride)
{
}
template <typename SplitPropertiesPortal>
VTKM_EXEC void operator()(const vtkm::Id& pointsToLeft,
const vtkm::Id& pointsToRight,
const vtkm::Range& lMaxRanges,
const vtkm::Range& rMinRanges,
const vtkm::FloatDefault& planeValue,
SplitPropertiesPortal& splits,
vtkm::Id inputIndex) const
{
SplitProperties split;
split.Plane = planeValue;
split.NumLeftPoints = pointsToLeft;
split.NumRightPoints = pointsToRight;
split.LMax = lMaxRanges.Max;
split.RMin = rMinRanges.Min;
split.Cost = vtkm::Abs(split.LMax * static_cast<vtkm::FloatDefault>(pointsToLeft) -
split.RMin * static_cast<vtkm::FloatDefault>(pointsToRight));
if (vtkm::IsNan(split.Cost))
{
split.Cost = vtkm::Infinity64();
}
splits.Set(inputIndex * Stride + Index, split);
//printf("Plane = %lf, NL = %lld, NR = %lld, LM = %lf, RM = %lf, C = %lf\n", split.Plane, split.NumLeftPoints, split.NumRightPoints, split.LMax, split.RMin, split.Cost);
}
vtkm::IdComponent Index;
vtkm::Id Stride;
};
struct SplitSelector : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>,
WholeArrayIn<>,
WholeArrayIn<>,
WholeArrayIn<>,
FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
using InputDomain = _1;
VTKM_CONT
SplitSelector(vtkm::IdComponent numPlanes,
vtkm::IdComponent maxLeafSize,
vtkm::IdComponent stride)
: NumPlanes(numPlanes)
, MaxLeafSize(maxLeafSize)
, Stride(stride)
{
}
template <typename SplitPropertiesPortal>
VTKM_EXEC void operator()(vtkm::Id index,
const SplitPropertiesPortal& xSplits,
const SplitPropertiesPortal& ySplits,
const SplitPropertiesPortal& zSplits,
const vtkm::Id& segmentSize,
TreeNode& node,
vtkm::FloatDefault& plane,
vtkm::Id& choice) const
{
if (segmentSize <= MaxLeafSize)
{
node.Dimension = -1;
choice = 0;
return;
}
choice = 1;
using Split = SplitProperties;
vtkm::FloatDefault minCost = vtkm::Infinity64();
const Split& xSplit = xSplits[ArgMin(xSplits, index * Stride, Stride)];
bool found = false;
if (xSplit.Cost < minCost && xSplit.NumLeftPoints != 0 && xSplit.NumRightPoints != 0)
{
minCost = xSplit.Cost;
node.Dimension = 0;
node.LMax = xSplit.LMax;
node.RMin = xSplit.RMin;
plane = xSplit.Plane;
found = true;
}
const Split& ySplit = ySplits[ArgMin(ySplits, index * Stride, Stride)];
if (ySplit.Cost < minCost && ySplit.NumLeftPoints != 0 && ySplit.NumRightPoints != 0)
{
minCost = ySplit.Cost;
node.Dimension = 1;
node.LMax = ySplit.LMax;
node.RMin = ySplit.RMin;
plane = ySplit.Plane;
found = true;
}
const Split& zSplit = zSplits[ArgMin(zSplits, index * Stride, Stride)];
if (zSplit.Cost < minCost && zSplit.NumLeftPoints != 0 && zSplit.NumRightPoints != 0)
{
minCost = zSplit.Cost;
node.Dimension = 2;
node.LMax = zSplit.LMax;
node.RMin = zSplit.RMin;
plane = zSplit.Plane;
found = true;
}
if (!found)
{
const Split& xMSplit = xSplits[NumPlanes];
minCost = xMSplit.Cost;
node.Dimension = 0;
node.LMax = xMSplit.LMax;
node.RMin = xMSplit.RMin;
plane = xMSplit.Plane;
const Split& yMSplit = ySplits[NumPlanes];
if (yMSplit.Cost < minCost && yMSplit.NumLeftPoints != 0 && yMSplit.NumRightPoints != 0)
{
minCost = yMSplit.Cost;
node.Dimension = 1;
node.LMax = yMSplit.LMax;
node.RMin = yMSplit.RMin;
plane = yMSplit.Plane;
}
const Split& zMSplit = zSplits[NumPlanes];
if (zMSplit.Cost < minCost && zMSplit.NumLeftPoints != 0 && zMSplit.NumRightPoints != 0)
{
minCost = zMSplit.Cost;
node.Dimension = 2;
node.LMax = zMSplit.LMax;
node.RMin = zMSplit.RMin;
plane = zMSplit.Plane;
}
}
//printf("Selected plane %lf, with cost %lf [%d, %lf, %lf]\n", plane, minCost, node.Dimension, node.LMax, node.RMin);
}
template <typename ArrayPortal>
VTKM_EXEC vtkm::Id ArgMin(const ArrayPortal& values, vtkm::Id start, vtkm::Id length) const
{
vtkm::Id minIdx = start;
for (vtkm::Id i = start; i < (start + length); ++i)
{
if (values[i].Cost < values[minIdx].Cost)
{
minIdx = i;
}
}
return minIdx;
}
vtkm::IdComponent NumPlanes;
vtkm::IdComponent MaxLeafSize;
vtkm::Id Stride;
};
struct CalculateSplitDirectionFlag : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::FloatDefault& x,
const vtkm::FloatDefault& y,
const vtkm::FloatDefault& z,
const TreeNode& split,
const vtkm::FloatDefault& plane,
vtkm::Id& flag) const
{
if (split.Dimension >= 0)
{
const vtkm::Vec<vtkm::FloatDefault, 3> point(x, y, z);
const vtkm::FloatDefault& c = point[split.Dimension];
// We use 0 to signify left child, 1 for right child
flag = 1 - static_cast<vtkm::Id>(c <= plane);
}
else
{
flag = 0;
}
}
}; // struct CalculateSplitDirectionFlag
struct SegmentSplitter : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4);
using InputDomain = _1;
VTKM_CONT
SegmentSplitter(vtkm::IdComponent maxLeafSize)
: MaxLeafSize(maxLeafSize)
{
}
VTKM_EXEC
void operator()(const vtkm::Id& segmentId,
const vtkm::Id& leqFlag,
const vtkm::Id& segmentSize,
vtkm::Id& newSegmentId) const
{
if (segmentSize <= MaxLeafSize)
{
// We do not split the segments which have cells fewer than MaxLeafSize, moving them to left
newSegmentId = 2 * segmentId;
}
else
{
newSegmentId = 2 * segmentId + leqFlag;
}
}
vtkm::IdComponent MaxLeafSize;
}; // struct SegmentSplitter
struct SplitIndicesCalculator : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::Id& leqFlag,
const vtkm::Id& trueFlagCount,
const vtkm::Id& countPreviousSegment,
const vtkm::Id& runningFalseFlagCount,
const vtkm::Id& totalFalseFlagCount,
vtkm::Id& scatterIndex) const
{
if (leqFlag)
{
scatterIndex = countPreviousSegment + totalFalseFlagCount + trueFlagCount;
}
else
{
scatterIndex = countPreviousSegment + runningFalseFlagCount - 1;
}
}
}; // struct SplitIndicesCalculator
struct Scatter : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>, FieldIn<>, WholeArrayOut<>);
typedef void ExecutionSignature(_1, _2, _3);
using InputDomain = _1;
template <typename InputType, typename OutputPortalType>
VTKM_EXEC void operator()(const InputType& in, const vtkm::Id& idx, OutputPortalType& out) const
{
out.Set(idx, in);
}
}; // struct Scatter
template <typename ValueArrayHandle, typename IndexArrayHandle>
ValueArrayHandle ScatterArray(const ValueArrayHandle& input, const IndexArrayHandle& indices)
{
ValueArrayHandle output;
output.Allocate(input.GetNumberOfValues());
vtkm::worklet::DispatcherMapField<Scatter>().Invoke(input, indices, output);
return output;
}
struct NonSplitIndexCalculator : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2);
using InputDomain = _1;
VTKM_CONT
NonSplitIndexCalculator(vtkm::IdComponent maxLeafSize)
: MaxLeafSize(maxLeafSize)
{
}
VTKM_EXEC void operator()(const vtkm::Id& inSegmentSize, vtkm::Id& outSegmentSize) const
{
if (inSegmentSize <= MaxLeafSize)
{
outSegmentSize = inSegmentSize;
}
else
{
outSegmentSize = 0;
}
}
vtkm::Id MaxLeafSize;
}; // struct NonSplitIndexCalculator
struct TreeLevelAdder : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldIn<>,
FieldIn<>,
FieldIn<>,
WholeArrayInOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
using InputDomain = _1;
VTKM_CONT
TreeLevelAdder(vtkm::Id cellIdsOffset, vtkm::Id treeOffset, vtkm::IdComponent maxLeafSize)
: CellIdsOffset(cellIdsOffset)
, TreeOffset(treeOffset)
, MaxLeafSize(maxLeafSize)
{
}
template <typename BoundingIntervalHierarchyPortal>
VTKM_EXEC void operator()(const vtkm::Id& index,
const TreeNode& split,
const vtkm::Id& start,
const vtkm::Id& count,
const vtkm::Id& numPreviousSplits,
BoundingIntervalHierarchyPortal& treePortal) const
{
vtkm:
cont::BoundingIntervalHierarchyNode node;
if (count > MaxLeafSize)
{