ConnectivityExplicitInternals.h 9.97 KB
Newer Older
1 2 3 4 5 6 7 8
//============================================================================
//  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.
//
9
//  Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 11 12
//  Copyright 2015 UT-Battelle, LLC.
//  Copyright 2015 Los Alamos National Security.
//
13
//  Under the terms of Contract DE-NA0003525 with NTESS,
14 15 16 17 18 19 20 21 22
//  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_cont_internal_ConnectivityExplicitInternals_h
#define vtk_m_cont_internal_ConnectivityExplicitInternals_h

23
#include <vtkm/CellShape.h>
24
#include <vtkm/cont/ArrayHandle.h>
25
#include <vtkm/cont/ArrayHandleCast.h>
26 27
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleCounting.h>
28
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
29
#include <vtkm/cont/internal/DeviceAdapterError.h>
30 31 32
#include <vtkm/exec/ExecutionWholeArray.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
33

34 35 36 37 38 39 40 41
namespace vtkm
{
namespace cont
{
namespace internal
{

template <typename NumIndicesArrayType, typename IndexOffsetArrayType, typename DeviceAdapterTag>
42 43 44 45
void buildIndexOffsets(const NumIndicesArrayType& numIndices,
                       IndexOffsetArrayType& offsets,
                       DeviceAdapterTag,
                       std::true_type)
46 47 48
{
  //We first need to make sure that NumIndices and IndexOffsetArrayType
  //have the same type so we can call scane exclusive
49
  using CastedNumIndicesType = vtkm::cont::ArrayHandleCast<vtkm::Id, NumIndicesArrayType>;
50 51 52 53

  // Although technically we are making changes to this object, the changes
  // are logically consistent with the previous state, so we consider it
  // valid under const.
54
  using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapterTag>;
55
  Algorithm::ScanExclusive(CastedNumIndicesType(numIndices), offsets);
56 57
}

58
template <typename NumIndicesArrayType, typename IndexOffsetArrayType, typename DeviceAdapterTag>
59 60 61
void buildIndexOffsets(const NumIndicesArrayType&,
                       IndexOffsetArrayType&,
                       DeviceAdapterTag,
62
                       std::false_type)
63 64 65 66 67 68 69 70
{
  //this is a no-op as the storage for the offsets is an implicit handle
  //and should already be built. This signature exists so that
  //the compiler doesn't try to generate un-used code that will
  //try and run Algorithm::ScanExclusive on an implicit array which will
  //cause a compile time failure.
}

71
template <typename ArrayHandleIndices, typename ArrayHandleOffsets, typename DeviceAdapterTag>
72 73
void buildIndexOffsets(const ArrayHandleIndices& numIndices,
                       ArrayHandleOffsets offsets,
74 75
                       DeviceAdapterTag tag)
{
76
  using IsWriteable = vtkm::cont::internal::IsWriteableArrayHandle<ArrayHandleOffsets>;
77 78 79
  buildIndexOffsets(numIndices, offsets, tag, typename IsWriteable::type());
}

80 81 82 83
template <typename ShapeStorageTag = VTKM_DEFAULT_STORAGE_TAG,
          typename NumIndicesStorageTag = VTKM_DEFAULT_STORAGE_TAG,
          typename ConnectivityStorageTag = VTKM_DEFAULT_STORAGE_TAG,
          typename IndexOffsetStorageTag = VTKM_DEFAULT_STORAGE_TAG>
84 85
struct ConnectivityExplicitInternals
{
86 87 88 89
  using ShapeArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8, ShapeStorageTag>;
  using NumIndicesArrayType = vtkm::cont::ArrayHandle<vtkm::IdComponent, NumIndicesStorageTag>;
  using ConnectivityArrayType = vtkm::cont::ArrayHandle<vtkm::Id, ConnectivityStorageTag>;
  using IndexOffsetArrayType = vtkm::cont::ArrayHandle<vtkm::Id, IndexOffsetStorageTag>;
90 91 92 93

  ShapeArrayType Shapes;
  NumIndicesArrayType NumIndices;
  ConnectivityArrayType Connectivity;
94
  mutable IndexOffsetArrayType IndexOffsets;
95 96

  bool ElementsValid;
97
  mutable bool IndexOffsetsValid;
98

99
  VTKM_CONT
100
  ConnectivityExplicitInternals()
101 102 103 104
    : ElementsValid(false)
    , IndexOffsetsValid(false)
  {
  }
105

106
  VTKM_CONT
107 108
  vtkm::Id GetNumberOfElements() const
  {
109 110
    VTKM_ASSERT(this->ElementsValid);

111 112 113
    return this->Shapes.GetNumberOfValues();
  }

114
  VTKM_CONT
115 116
  void ReleaseResourcesExecution()
  {
117 118 119 120 121 122
    this->Shapes.ReleaseResourcesExecution();
    this->NumIndices.ReleaseResourcesExecution();
    this->Connectivity.ReleaseResourcesExecution();
    this->IndexOffsets.ReleaseResourcesExecution();
  }

123 124
  template <typename Device>
  VTKM_CONT void BuildIndexOffsets(Device) const
125
  {
126
    VTKM_ASSERT(this->ElementsValid);
127

128
    if (!this->IndexOffsetsValid)
129
    {
130
      buildIndexOffsets(this->NumIndices, this->IndexOffsets, Device());
131
      this->IndexOffsetsValid = true;
132
    }
133 134
  }

135
  VTKM_CONT
136 137 138 139
  void BuildIndexOffsets(vtkm::cont::DeviceAdapterTagError) const
  {
    if (!this->IndexOffsetsValid)
    {
140
      throw vtkm::cont::ErrorBadType(
141
        "Cannot build indices using the error device. Must be created previously.");
142 143 144
    }
  }

145
  VTKM_CONT
146
  void PrintSummary(std::ostream& out) const
147
  {
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
    if (this->ElementsValid)
    {
      out << "     Shapes: ";
      vtkm::cont::printSummary_ArrayHandle(this->Shapes, out);
      out << "     NumIndices: ";
      vtkm::cont::printSummary_ArrayHandle(this->NumIndices, out);
      out << "     Connectivity: ";
      vtkm::cont::printSummary_ArrayHandle(this->Connectivity, out);
      if (this->IndexOffsetsValid)
      {
        out << "     IndexOffsets: ";
        vtkm::cont::printSummary_ArrayHandle(this->IndexOffsets, out);
      }
      else
      {
        out << "     IndexOffsets: Not Allocated" << std::endl;
      }
    }
    else
    {
      out << "     Not Allocated" << std::endl;
    }
170 171
  }
};
172 173 174 175 176 177


// Worklet to expand the PointToCell numIndices array by repeating cell index
class ExpandIndices : public vtkm::worklet::WorkletMapField
{
public:
178
  using ControlSignature = void(FieldIn<> cellIndex,
179 180 181
                                FieldIn<> offset,
                                FieldIn<> numIndices,
                                WholeArrayOut<> cellIndices);
182
  using ExecutionSignature = void(_1, _2, _3, _4);
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
  using InputDomain = _1;

  VTKM_CONT
  ExpandIndices() {}

  template <typename PortalType>
  VTKM_EXEC void operator()(const vtkm::Id& cellIndex,
                            const vtkm::Id& offset,
                            const vtkm::Id& numIndices,
                            const PortalType& cellIndices) const
  {
    VTKM_ASSERT(cellIndices.GetNumberOfValues() >= offset + numIndices);
    vtkm::Id startIndex = offset;
    for (vtkm::Id i = 0; i < numIndices; i++)
    {
      cellIndices.Set(startIndex++, cellIndex);
    }
  }
};

class ScatterValues : public vtkm::worklet::WorkletMapField
{
public:
206 207
  using ControlSignature = void(FieldIn<> index, FieldIn<> value, WholeArrayOut<> output);
  using ExecutionSignature = void(_1, _2, _3);
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
  using InputDomain = _1;

  template <typename T, typename PortalType>
  VTKM_EXEC void operator()(const vtkm::Id& index, const T& value, const PortalType& output) const
  {
    output.Set(index, value);
  }
};

template <typename PointToCell, typename C2PShapeStorageTag, typename Device>
void ComputeCellToPointConnectivity(ConnectivityExplicitInternals<C2PShapeStorageTag>& cell2Point,
                                    const PointToCell& point2Cell,
                                    vtkm::Id numberOfPoints,
                                    Device)
{
  // PointToCell connectivity array (point indices) will be
  // transformed into the CellToPoint numIndices array using reduction
  //
  // PointToCell numIndices array using expansion will be
  // transformed into the CellToPoint connectivity array

  if (cell2Point.ElementsValid)
  {
    return;
  }

  using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<Device>;

  // Sizes of the PointToCell information
  vtkm::Id numberOfCells = point2Cell.NumIndices.GetNumberOfValues();
  vtkm::Id connectivityLength = point2Cell.Connectivity.GetNumberOfValues();

  // PointToCell connectivity will be basis of CellToPoint numIndices
  vtkm::cont::ArrayHandle<vtkm::Id> pointIndices;
  Algorithm::Copy(point2Cell.Connectivity, pointIndices);

  // PointToCell numIndices will be basis of CellToPoint connectivity

  cell2Point.Connectivity.Allocate(connectivityLength);
  vtkm::cont::ArrayHandleCounting<vtkm::Id> index(0, 1, numberOfCells);

  vtkm::worklet::DispatcherMapField<ExpandIndices, Device> expandDispatcher;
  expandDispatcher.Invoke(
    index, point2Cell.IndexOffsets, point2Cell.NumIndices, cell2Point.Connectivity);

  // SortByKey where key is PointToCell connectivity and value is the expanded cellIndex
  Algorithm::SortByKey(pointIndices, cell2Point.Connectivity);

  // CellToPoint numIndices from the now sorted PointToCell connectivity
  vtkm::cont::ArrayHandleConstant<vtkm::IdComponent> numArray(1, connectivityLength);
  vtkm::cont::ArrayHandle<vtkm::Id> uniquePoints;
  vtkm::cont::ArrayHandle<vtkm::IdComponent> numIndices;
  Algorithm::ReduceByKey(pointIndices, numArray, uniquePoints, numIndices, vtkm::Add());

  // if not all the points have a cell
  if (uniquePoints.GetNumberOfValues() < numberOfPoints)
  {
    vtkm::cont::ArrayHandle<vtkm::IdComponent> fullNumIndices;
    Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::IdComponent>(0, numberOfPoints),
                    fullNumIndices);
    vtkm::worklet::DispatcherMapField<ScatterValues, Device>().Invoke(
      uniquePoints, numIndices, fullNumIndices);
    numIndices = fullNumIndices;
  }

  // Set the CellToPoint information
  cell2Point.Shapes = vtkm::cont::make_ArrayHandleConstant(
    static_cast<vtkm::UInt8>(CELL_SHAPE_VERTEX), numberOfPoints);
  cell2Point.NumIndices = numIndices;

  cell2Point.ElementsValid = true;
  cell2Point.IndexOffsetsValid = false;
}
281 282 283 284 285
}
}
} // namespace vtkm::cont::internal

#endif //vtk_m_cont_internal_ConnectivityExplicitInternals_h