Add OrientNormals worklet.

vtkm/worklet/CMakeLists.txt

@@ -49,6 +49,10 @@ set(headers
   NDimsHistMarginalization.h
   NDimsHistogram.h
   Normalize.h
+  OrientCellNormals.h
+  OrientNormals.h
+  OrientPointNormals.h
+  OrientPointAndCellNormals.h
   OscillatorSource.h
   ParticleAdvection.h
   PointAverage.h

vtkm/worklet/OrientCellNormals.h

@@ -0,0 +1,451 @@
+//============================================================================
+//  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 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2019 UT-Battelle, LLC.
+//  Copyright 2019 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 vtkm_m_worklet_OrientCellNormals_h
+#define vtkm_m_worklet_OrientCellNormals_h
+
+#include <vtkm/Range.h>
+#include <vtkm/Types.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/ArrayHandleBitField.h>
+#include <vtkm/cont/ArrayHandleConstant.h>
+#include <vtkm/cont/ArrayHandleTransform.h>
+#include <vtkm/cont/ArrayRangeCompute.h>
+#include <vtkm/cont/BitField.h>
+#include <vtkm/cont/Logging.h>
+
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/MaskIndices.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+
+#include <vtkm/VecTraits.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+
+///
+/// Orients normals to point outside of the dataset. This requires a closed
+/// manifold surface or else the behavior is undefined. This requires an
+/// unstructured cellset as input.
+///
+class OrientCellNormals
+{
+  static constexpr vtkm::Id INVALID_ID = -1;
+
+  // Returns true if v1 and v2 are pointing in the same hemisphere.
+  template <typename T>
+  VTKM_EXEC static bool SameDirection(const vtkm::Vec<T, 3>& v1, const vtkm::Vec<T, 3>& v2)
+  {
+    return vtkm::Dot(v1, v2) >= 0;
+  }
+
+  // Ensure that the normal is pointing in the same hemisphere as ref.
+  // Returns true if normal is modified.
+  template <typename T>
+  VTKM_EXEC static bool Align(vtkm::Vec<T, 3>& normal, const vtkm::Vec<T, 3>& ref)
+  {
+    if (!SameDirection(normal, ref))
+    {
+      normal = -normal;
+      return true;
+    }
+    return false;
+  }
+
+public:
+  // Locates starting points for BFS traversal of dataset by finding points
+  // on the dataset boundaries. These points are marked as active.
+  // Initializes the ActivePoints array.
+  class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldIn coords, WholeArrayIn ranges, FieldOut activePoints);
+    using ExecutionSignature = _3(_1 coord, _2 ranges);
+
+    template <typename CoordT, typename RangePortal>
+    VTKM_EXEC bool operator()(const vtkm::Vec<CoordT, 3>& point, const RangePortal& ranges) const
+    {
+      bool isActive = false;
+      for (vtkm::IdComponent dim = 0; dim < 3; ++dim)
+      {
+        const auto& range = ranges.Get(dim);
+        const auto val = static_cast<decltype(range.Min)>(point[dim]);
+        if (val <= range.Min || val >= range.Max)
+        {
+          isActive = true;
+        }
+      }
+      return isActive;
+    }
+  };
+
+  // For each of the source points, determine the boundaries it lies on. Align
+  // each incident cell's normal to point out of the boundary, marking each cell
+  // as both visited and active.
+  // Clears the active flags for points, and marks the current point as visited.
+  class WorkletProcessSourceCells : public vtkm::worklet::WorkletMapCellToPoint
+  {
+  public:
+    using ControlSignature = void(CellSetIn cells,
+                                  FieldInPoint coords,
+                                  WholeArrayIn ranges,
+                                  WholeArrayInOut cellNormals,
+                                  // InOut for preserve data on masked indices
+                                  BitFieldInOut activeCells,
+                                  BitFieldInOut visitedCells,
+                                  FieldInOutPoint activePoints,
+                                  FieldInOutPoint visitedPoints);
+    using ExecutionSignature = void(CellIndices cellIds,
+                                    _2 coords,
+                                    _3 ranges,
+                                    _4 cellNormals,
+                                    _5 activeCells,
+                                    _6 visitedCells,
+                                    _7 activePoints,
+                                    _8 visitedPoints);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename CellList,
+              typename CoordComp,
+              typename RangePortal,
+              typename CellNormalPortal,
+              typename ActiveCellsBitPortal,
+              typename VisitedCellsBitPortal>
+    VTKM_EXEC void operator()(const CellList& cellIds,
+                              const vtkm::Vec<CoordComp, 3>& coord,
+                              const RangePortal& ranges,
+                              CellNormalPortal& cellNormals,
+                              ActiveCellsBitPortal& activeCells,
+                              VisitedCellsBitPortal& visitedCells,
+                              bool& pointIsActive,
+                              bool& pointIsVisited) const
+    {
+      using NormalType = typename CellNormalPortal::ValueType;
+      VTKM_STATIC_ASSERT_MSG(vtkm::VecTraits<NormalType>::NUM_COMPONENTS == 3,
+                             "Cell Normals expected to have 3 components.");
+      using NormalCompType = typename NormalType::ComponentType;
+
+      // Find the vector that points out of the dataset from the current point:
+      const NormalType refNormal = [&]() -> NormalType {
+        NormalType normal{ NormalCompType{ 0 } };
+        NormalCompType numNormals{ 0 };
+        for (vtkm::IdComponent dim = 0; dim < 3; ++dim)
+        {
+          const auto range = ranges.Get(dim);
+          const auto val = static_cast<decltype(range.Min)>(coord[dim]);
+          if (val <= range.Min)
+          {
+            NormalType compNormal{ NormalCompType{ 0 } };
+            compNormal[dim] = NormalCompType{ -1 };
+            normal += compNormal;
+            ++numNormals;
+          }
+          else if (val >= range.Max)
+          {
+            NormalType compNormal{ NormalCompType{ 0 } };
+            compNormal[dim] = NormalCompType{ 1 };
+            normal += compNormal;
+            ++numNormals;
+          }
+        }
+
+        VTKM_ASSERT("Source point is not on a boundary?" && numNormals > 0.5);
+        return normal / numNormals;
+      }();
+
+      // Align all cell normals to the reference, marking them as active and
+      // visited.
+      const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents();
+      for (vtkm::IdComponent c = 0; c < numCells; ++c)
+      {
+        const vtkm::Id cellId = cellIds[c];
+
+        if (!visitedCells.OrBitAtomic(cellId, true))
+        { // This thread is the first to touch this cell.
+          activeCells.SetBitAtomic(cellId, true);
+
+          NormalType cellNormal = cellNormals.Get(cellId);
+          if (Align(cellNormal, refNormal))
+          {
+            cellNormals.Set(cellId, cellNormal);
+          }
+        }
+      }
+
+      // Mark current point as inactive but visited:
+      pointIsActive = false;
+      pointIsVisited = true;
+    }
+  };
+
+  // Mark each incident point as active and visited.
+  // Marks the current cell as inactive.
+  class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell
+  {
+  public:
+    using ControlSignature = void(CellSetIn cell,
+                                  BitFieldInOut activePoints,
+                                  BitFieldInOut visitedPoints,
+                                  FieldInOutCell activeCells);
+    using ExecutionSignature = _4(PointIndices pointIds, _2 activePoints, _3 visitedPoints);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename PointList, typename ActivePointsBitPortal, typename VisitedPointsBitPortal>
+    VTKM_EXEC bool operator()(const PointList& pointIds,
+                              ActivePointsBitPortal& activePoints,
+                              VisitedPointsBitPortal& visitedPoints) const
+    {
+      const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents();
+      for (vtkm::IdComponent p = 0; p < numPoints; ++p)
+      {
+        const vtkm::Id pointId = pointIds[p];
+        if (!visitedPoints.OrBitAtomic(pointId, true))
+        { // This thread owns this point.
+          activePoints.SetBitAtomic(pointId, true);
+        }
+      }
+
+      // Mark current cell as inactive:
+      return false;
+    }
+  };
+
+  // Mark each incident cell as active, setting a visited neighbor
+  // cell as its reference for alignment.
+  // Marks the current point as inactive.
+  class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint
+  {
+  public:
+    using ControlSignature = void(CellSetIn cells,
+                                  WholeArrayOut refCells,
+                                  BitFieldInOut activeCells,
+                                  BitFieldIn visitedCells,
+                                  FieldInOutPoint activePoints);
+    using ExecutionSignature = _5(CellIndices cellIds,
+                                  _2 refCells,
+                                  _3 activeCells,
+                                  _4 visitedCells);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename CellList,
+              typename RefCellPortal,
+              typename ActiveCellBitPortal,
+              typename VisitedCellBitPortal>
+    VTKM_EXEC bool operator()(const CellList& cellIds,
+                              RefCellPortal& refCells,
+                              ActiveCellBitPortal& activeCells,
+                              const VisitedCellBitPortal& visitedCells) const
+    {
+      // One of the cells must be marked visited already. Find it and use it as
+      // an alignment reference for the others:
+      const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents();
+      const vtkm::Id refCellId = [&]() -> vtkm::Id {
+        for (vtkm::IdComponent c = 0; c < numCells; ++c)
+        {
+          const vtkm::Id cellId = cellIds[c];
+          if (visitedCells.GetBit(cellId))
+          {
+            return cellId;
+          }
+        }
+        return INVALID_ID;
+      }();
+
+      VTKM_ASSERT("No reference cell found." && refCellId != INVALID_ID);
+
+      for (vtkm::IdComponent c = 0; c < numCells; ++c)
+      {
+        const vtkm::Id cellId = cellIds[c];
+        if (!visitedCells.GetBit(cellId))
+        {
+          if (!activeCells.OrBitAtomic(cellId, true))
+          { // This thread owns this cell.
+            refCells.Set(cellId, refCellId);
+          }
+        }
+      }
+
+      // Mark current point as inactive:
+      return false;
+    }
+  };
+
+  // Align the normal of each active cell, to its reference cell normal.
+  // The cell is marked visited.
+  class WorkletProcessCellNormals : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldIn refCells,
+                                  WholeArrayInOut cellNormals,
+                                  FieldInOut visitedCells);
+    using ExecutionSignature = _3(InputIndex cellId, _1 refCellId, _2 cellNormals);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename CellNormalsPortal>
+    VTKM_EXEC bool operator()(const vtkm::Id cellId,
+                              const vtkm::Id refCellId,
+                              CellNormalsPortal& cellNormals) const
+    {
+      const auto refNormal = cellNormals.Get(refCellId);
+      auto normal = cellNormals.Get(cellId);
+      if (Align(normal, refNormal))
+      {
+        cellNormals.Set(cellId, normal);
+      }
+
+      // Mark cell as visited:
+      return true;
+    }
+  };
+
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename CellNormalCompType,
+            typename CellNormalStorageType>
+  VTKM_CONT static void Run(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<CellNormalCompType, 3>, CellNormalStorageType>& cellNormals)
+  {
+    using RangeType = vtkm::cont::ArrayHandle<vtkm::Range>;
+
+    using MarkSourcePoints = vtkm::worklet::DispatcherMapField<WorkletMarkSourcePoints>;
+    using ProcessSourceCells = vtkm::worklet::DispatcherMapTopology<WorkletProcessSourceCells>;
+    using MarkActivePoints = vtkm::worklet::DispatcherMapTopology<WorkletMarkActivePoints>;
+    using MarkActiveCells = vtkm::worklet::DispatcherMapTopology<WorkletMarkActiveCells>;
+    using ProcessCellNormals = vtkm::worklet::DispatcherMapField<WorkletProcessCellNormals>;
+
+    const vtkm::Id numPoints = coords.GetNumberOfValues();
+    const vtkm::Id numCells = cells.GetNumberOfCells();
+
+    VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf,
+                   "OrientCellNormals worklet (%lld points, %lld cells)",
+                   static_cast<vtkm::Int64>(numPoints),
+                   static_cast<vtkm::Int64>(numCells));
+
+    // active = cells / point to be used in the next worklet invocation mask.
+    vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints
+    auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits);
+
+    vtkm::cont::BitField activeCellBits;
+    vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells);
+    auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits);
+
+    // visited = cells / points that have been corrected.
+    vtkm::cont::BitField visitedPointBits;
+    vtkm::cont::Algorithm::Fill(visitedPointBits, false, numPoints);
+    auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits);
+
+    vtkm::cont::BitField visitedCellBits;
+    vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells);
+    auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits);
+
+    vtkm::cont::ArrayHandle<vtkm::Id> mask; // Allocated as needed
+
+    // For each cell, store a reference alignment cell.
+    vtkm::cont::ArrayHandle<vtkm::Id> refCells;
+    {
+      vtkm::cont::Algorithm::Copy(
+        vtkm::cont::make_ArrayHandleConstant<vtkm::Id>(INVALID_ID, numCells), refCells);
+    }
+
+    // 1) Compute range of coords.
+    const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords);
+
+    // 2) Locate points on a boundary, since their normal alignment direction
+    //    is known.
+    {
+      MarkSourcePoints dispatcher;
+      dispatcher.Invoke(coords, ranges, activePoints);
+    }
+
+    // 3) For each source point, align the normals of the adjacent cells.
+    {
+      vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+      (void)numActive;
+      VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                 "ProcessSourceCells from " << numActive << " source points.");
+      ProcessSourceCells dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+      dispatcher.Invoke(cells,
+                        coords,
+                        ranges,
+                        cellNormals,
+                        activeCellBits,
+                        visitedCellBits,
+                        activePoints,
+                        visitedPoints);
+    }
+
+    for (size_t iter = 1;; ++iter)
+    {
+      // 4) Mark unvisited points adjacent to active cells
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActivePoints from " << numActive << " active cells.");
+        MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, activePointBits, visitedPointBits, activeCells);
+      }
+
+      // 5) Mark unvisited cells adjacent to active points
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActiveCells from " << numActive << " active points.");
+        MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, refCells, activeCellBits, visitedCellBits, activePoints);
+      }
+
+      vtkm::Id numActiveCells = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask);
+
+      if (numActiveCells == 0)
+      { // Done!
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf, "Iteration " << iter << ": Traversal complete.");
+        break;
+      }
+
+      VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                 "Iteration " << iter << ": Processing " << numActiveCells << " normals.");
+
+      // 5) Correct normals for active cells.
+      {
+        ProcessCellNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(refCells, cellNormals, visitedCells);
+      }
+    }
+  }
+};
+}
+} // end namespace vtkm::worklet
+
+
+#endif // vtkm_m_worklet_OrientCellNormals_h

vtkm/worklet/OrientNormals.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 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2019 UT-Battelle, LLC.
+//  Copyright 2019 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 vtkm_m_worklet_OrientNormals_h
+#define vtkm_m_worklet_OrientNormals_h
+
+#include <vtkm/Types.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/ArrayHandleTransform.h>
+
+#include <vtkm/worklet/OrientCellNormals.h>
+#include <vtkm/worklet/OrientPointAndCellNormals.h>
+#include <vtkm/worklet/OrientPointNormals.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+
+///
+/// Orients normals to point outside of the dataset. This requires a closed
+/// manifold surface or else the behavior is undefined. This requires an
+/// unstructured cellset as input.
+///
+class OrientNormals
+{
+public:
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename CellNormalCompType,
+            typename CellNormalStorageType>
+  VTKM_CONT static void RunCellNormals(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<CellNormalCompType, 3>, CellNormalStorageType>& cellNormals)
+  {
+    OrientCellNormals::Run(cells, coords, cellNormals);
+  }
+
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename PointNormalCompType,
+            typename PointNormalStorageType>
+  VTKM_CONT static void RunPointNormals(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<PointNormalCompType, 3>, PointNormalStorageType>&
+      pointNormals)
+  {
+    OrientPointNormals::Run(cells, coords, pointNormals);
+  }
+
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename PointNormalCompType,
+            typename PointNormalStorageType,
+            typename CellNormalCompType,
+            typename CellNormalStorageType>
+  VTKM_CONT static void RunPointAndCellNormals(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<PointNormalCompType, 3>, PointNormalStorageType>&
+      pointNormals,
+    vtkm::cont::ArrayHandle<vtkm::Vec<CellNormalCompType, 3>, CellNormalStorageType>& cellNormals)
+  {
+    OrientPointAndCellNormals::Run(cells, coords, pointNormals, cellNormals);
+  }
+
+  struct NegateFunctor
+  {
+    template <typename T>
+    VTKM_EXEC_CONT T operator()(const T& val) const
+    {
+      return -val;
+    }
+  };
+
+  ///
+  /// Reverse the normals to point in the opposite direction.
+  ///
+  template <typename NormalCompType, typename NormalStorageType>
+  VTKM_CONT static void RunFlipNormals(
+    vtkm::cont::ArrayHandle<vtkm::Vec<NormalCompType, 3>, NormalStorageType>& normals)
+  {
+    const auto flippedAlias = vtkm::cont::make_ArrayHandleTransform(normals, NegateFunctor{});
+    vtkm::cont::Algorithm::Copy(flippedAlias, normals);
+  }
+};
+}
+} // end namespace vtkm::worklet
+
+
+#endif // vtkm_m_worklet_OrientNormals_h

vtkm/worklet/OrientPointAndCellNormals.h

@@ -0,0 +1,421 @@
+//============================================================================
+//  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 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2019 UT-Battelle, LLC.
+//  Copyright 2019 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 vtkm_m_worklet_OrientPointAndCellNormals_h
+#define vtkm_m_worklet_OrientPointAndCellNormals_h
+
+#include <vtkm/Range.h>
+#include <vtkm/Types.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/ArrayHandleBitField.h>
+#include <vtkm/cont/ArrayHandleConstant.h>
+#include <vtkm/cont/ArrayHandleTransform.h>
+#include <vtkm/cont/ArrayRangeCompute.h>
+#include <vtkm/cont/BitField.h>
+#include <vtkm/cont/Logging.h>
+
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/MaskIndices.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+
+#include <vtkm/VecTraits.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+
+///
+/// Orients normals to point outside of the dataset. This requires a closed
+/// manifold surface or else the behavior is undefined. This requires an
+/// unstructured cellset as input.
+///
+class OrientPointAndCellNormals
+{
+  static constexpr vtkm::Id INVALID_ID = -1;
+
+  // Returns true if v1 and v2 are pointing in the same hemisphere.
+  template <typename T>
+  VTKM_EXEC static bool SameDirection(const vtkm::Vec<T, 3>& v1, const vtkm::Vec<T, 3>& v2)
+  {
+    return vtkm::Dot(v1, v2) >= 0;
+  }
+
+  // Ensure that the normal is pointing in the same hemisphere as ref.
+  // Returns true if normal is modified.
+  template <typename T>
+  VTKM_EXEC static bool Align(vtkm::Vec<T, 3>& normal, const vtkm::Vec<T, 3>& ref)
+  {
+    if (!SameDirection(normal, ref))
+    {
+      normal = -normal;
+      return true;
+    }
+    return false;
+  }
+
+public:
+  // Locates starting points for BFS traversal of dataset by finding points
+  // on the dataset boundaries. The normals for these points are corrected by
+  // making them point outside of the dataset, and they are marked as both
+  // active and visited.
+  class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldIn coords,
+                                  FieldInOut pointNormals,
+                                  WholeArrayIn ranges,
+                                  FieldOut activePoints,
+                                  FieldOut visitedPoints);
+    using ExecutionSignature =
+      void(_1 coord, _2 pointNormal, _3 ranges, _4 activePoints, _5 visitedPoints);
+
+    template <typename CoordT, typename NormalT, typename RangePortal>
+    VTKM_EXEC void operator()(const vtkm::Vec<CoordT, 3>& point,
+                              vtkm::Vec<NormalT, 3>& pointNormal,
+                              const RangePortal& ranges,
+                              bool& isActive,
+                              bool& isVisited) const
+    {
+      for (vtkm::IdComponent dim = 0; dim < 3; ++dim)
+      {
+        const auto& range = ranges.Get(dim);
+        const auto val = static_cast<decltype(range.Min)>(point[dim]);
+        if (val <= range.Min)
+        {
+          vtkm::Vec<NormalT, 3> ref{ static_cast<NormalT>(0) };
+          ref[dim] = static_cast<NormalT>(-1);
+          Align(pointNormal, ref);
+          isActive = true;
+          isVisited = true;
+          return;
+        }
+        else if (val >= range.Max)
+        {
+          vtkm::Vec<NormalT, 3> ref{ static_cast<NormalT>(0) };
+          ref[dim] = static_cast<NormalT>(1);
+          Align(pointNormal, ref);
+          isActive = true;
+          isVisited = true;
+          return;
+        }
+      }
+
+      isActive = false;
+      isVisited = false;
+    }
+  };
+
+  // Mark each incident cell as active and visited.
+  // Marks the current point as inactive.
+  class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint
+  {
+  public:
+    using ControlSignature = void(CellSetIn cell,
+                                  BitFieldInOut activeCells,
+                                  BitFieldInOut visitedCells,
+                                  FieldInOutPoint activePoints);
+    using ExecutionSignature = _4(CellIndices cellIds, _2 activeCells, _3 visitedCells);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename CellList, typename ActiveCellsBitPortal, typename VisitedCellsBitPortal>
+    VTKM_EXEC bool operator()(const CellList& cellIds,
+                              ActiveCellsBitPortal& activeCells,
+                              VisitedCellsBitPortal& visitedCells) const
+    {
+      const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents();
+      for (vtkm::IdComponent c = 0; c < numCells; ++c)
+      {
+        const vtkm::Id cellId = cellIds[c];
+        if (!visitedCells.OrBitAtomic(cellId, true))
+        { // This thread owns this cell.
+          activeCells.SetBitAtomic(cellId, true);
+        }
+      }
+
+      // Mark current point as inactive:
+      return false;
+    }
+  };
+
+  // Align the current cell's normals to an adjacent visited point's normal.
+  class WorkletProcessCellNormals : public vtkm::worklet::WorkletMapPointToCell
+  {
+  public:
+    using ControlSignature = void(CellSetIn cells,
+                                  WholeArrayIn pointNormals,
+                                  WholeArrayInOut cellNormals,
+                                  BitFieldIn visitedPoints);
+    using ExecutionSignature = void(PointIndices pointIds,
+                                    InputIndex cellId,
+                                    _2 pointNormals,
+                                    _3 cellNormals,
+                                    _4 visitedPoints);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename PointList,
+              typename PointNormalsPortal,
+              typename CellNormalsPortal,
+              typename VisitedPointsBitPortal>
+    VTKM_EXEC void operator()(const PointList& pointIds,
+                              const vtkm::Id cellId,
+                              const PointNormalsPortal& pointNormals,
+                              CellNormalsPortal& cellNormals,
+                              const VisitedPointsBitPortal& visitedPoints) const
+    {
+      // Use the normal of a visited point as a reference:
+      const vtkm::Id refPointId = [&]() -> vtkm::Id {
+        const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents();
+        for (vtkm::IdComponent p = 0; p < numPoints; ++p)
+        {
+          const vtkm::Id pointId = pointIds[p];
+          if (visitedPoints.GetBit(pointId))
+          {
+            return pointId;
+          }
+        }
+
+        return INVALID_ID;
+      }();
+
+      VTKM_ASSERT("No reference point found." && refPointId != INVALID_ID);
+
+      const auto refNormal = pointNormals.Get(refPointId);
+      auto normal = cellNormals.Get(cellId);
+      if (Align(normal, refNormal))
+      {
+        cellNormals.Set(cellId, normal);
+      }
+    }
+  };
+
+  // Mark each incident point as active and visited.
+  // Marks the current cell as inactive.
+  class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell
+  {
+  public:
+    using ControlSignature = void(CellSetIn cell,
+                                  BitFieldInOut activePoints,
+                                  BitFieldInOut visitedPoints,
+                                  FieldInOutCell activeCells);
+    using ExecutionSignature = _4(PointIndices pointIds, _2 activePoint, _3 visitedPoint);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename PointList, typename ActivePointsBitPortal, typename VisitedPointsBitPortal>
+    VTKM_EXEC bool operator()(const PointList& pointIds,
+                              ActivePointsBitPortal& activePoints,
+                              VisitedPointsBitPortal& visitedPoints) const
+    {
+      const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents();
+      for (vtkm::IdComponent p = 0; p < numPoints; ++p)
+      {
+        const vtkm::Id pointId = pointIds[p];
+        if (!visitedPoints.OrBitAtomic(pointId, true))
+        { // This thread owns this point.
+          activePoints.SetBitAtomic(pointId, true);
+        }
+      }
+
+      // Mark current cell as inactive:
+      return false;
+    }
+  };
+
+  // Align the current point's normals to an adjacent visited cell's normal.
+  class WorkletProcessPointNormals : public vtkm::worklet::WorkletMapCellToPoint
+  {
+  public:
+    using ControlSignature = void(CellSetIn cells,
+                                  WholeArrayInOut pointNormals,
+                                  WholeArrayIn cellNormals,
+                                  BitFieldIn visitedCells);
+    using ExecutionSignature = void(CellIndices cellIds,
+                                    InputIndex pointId,
+                                    _2 pointNormals,
+                                    _3 cellNormals,
+                                    _4 visitedCells);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename CellList,
+              typename CellNormalsPortal,
+              typename PointNormalsPortal,
+              typename VisitedCellsBitPortal>
+    VTKM_EXEC void operator()(const CellList& cellIds,
+                              const vtkm::Id pointId,
+                              PointNormalsPortal& pointNormals,
+                              const CellNormalsPortal& cellNormals,
+                              const VisitedCellsBitPortal& visitedCells) const
+    {
+      // Use the normal of a visited cell as a reference:
+      const vtkm::Id refCellId = [&]() -> vtkm::Id {
+        const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents();
+        for (vtkm::IdComponent c = 0; c < numCells; ++c)
+        {
+          const vtkm::Id cellId = cellIds[c];
+          if (visitedCells.GetBit(cellId))
+          {
+            return cellId;
+          }
+        }
+
+        return INVALID_ID;
+      }();
+
+      VTKM_ASSERT("No reference cell found." && refCellId != INVALID_ID);
+
+      const auto refNormal = cellNormals.Get(refCellId);
+      auto normal = pointNormals.Get(pointId);
+      if (Align(normal, refNormal))
+      {
+        pointNormals.Set(pointId, normal);
+      }
+    }
+  };
+
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename PointNormalCompType,
+            typename PointNormalStorageType,
+            typename CellNormalCompType,
+            typename CellNormalStorageType>
+  VTKM_CONT static void Run(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<PointNormalCompType, 3>, PointNormalStorageType>&
+      pointNormals,
+    vtkm::cont::ArrayHandle<vtkm::Vec<CellNormalCompType, 3>, CellNormalStorageType>& cellNormals)
+  {
+    using RangeType = vtkm::cont::ArrayHandle<vtkm::Range>;
+
+    using MarkSourcePoints = vtkm::worklet::DispatcherMapField<WorkletMarkSourcePoints>;
+    using MarkActiveCells = vtkm::worklet::DispatcherMapTopology<WorkletMarkActiveCells>;
+    using ProcessCellNormals = vtkm::worklet::DispatcherMapTopology<WorkletProcessCellNormals>;
+    using MarkActivePoints = vtkm::worklet::DispatcherMapTopology<WorkletMarkActivePoints>;
+    using ProcessPointNormals = vtkm::worklet::DispatcherMapTopology<WorkletProcessPointNormals>;
+
+    const vtkm::Id numCells = cells.GetNumberOfCells();
+
+    VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf,
+                   "OrientPointAndCellNormals worklet (%lld points, %lld cells)",
+                   static_cast<vtkm::Int64>(coords.GetNumberOfValues()),
+                   static_cast<vtkm::Int64>(numCells));
+
+    // active = cells / point to be used in the next worklet invocation mask.
+    vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints
+    auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits);
+
+    vtkm::cont::BitField activeCellBits;
+    vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells);
+    auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits);
+
+    // visited = cells / points that have been corrected.
+    vtkm::cont::BitField visitedPointBits; // Initialized by MarkSourcePoints
+    auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits);
+
+    vtkm::cont::BitField visitedCellBits;
+    vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells);
+    auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits);
+
+    vtkm::cont::ArrayHandle<vtkm::Id> mask; // Allocated as needed
+
+    // 1) Compute range of coords.
+    const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords);
+
+    // 2) Locate points on a boundary and align their normal to point out of the
+    //    dataset:
+    {
+      MarkSourcePoints dispatcher;
+      dispatcher.Invoke(coords, pointNormals, ranges, activePoints, visitedPoints);
+    }
+
+    for (size_t iter = 1;; ++iter)
+    {
+      // 3) Mark unvisited cells adjacent to active points
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActiveCells from " << numActive << " active points.");
+        MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, activeCellBits, visitedCellBits, activePoints);
+      }
+
+      vtkm::Id numActiveCells = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask);
+
+      if (numActiveCells == 0)
+      { // Done!
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "Iteration " << iter << ": Traversal complete; no more cells");
+        break;
+      }
+
+      VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                 "Iteration " << iter << ": Processing " << numActiveCells << " cell normals.");
+
+      // 4) Correct normals for active cells.
+      {
+        ProcessCellNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, pointNormals, cellNormals, visitedPointBits);
+      }
+
+      // 5) Mark unvisited points adjacent to active cells
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActivePoints from " << numActive << " active cells.");
+        MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, activePointBits, visitedPointBits, activeCells);
+      }
+
+      vtkm::Id numActivePoints =
+        vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+
+      if (numActivePoints == 0)
+      { // Done!
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "Iteration " << iter << ": Traversal complete; no more points");
+        break;
+      }
+
+      VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                 "Iteration " << iter << ": Processing " << numActivePoints << " point normals.");
+
+      // 4) Correct normals for active points.
+      {
+        ProcessPointNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, pointNormals, cellNormals, visitedCellBits);
+      }
+    }
+  }
+};
+}
+} // end namespace vtkm::worklet
+
+
+#endif // vtkm_m_worklet_OrientPointAndCellNormals_h

vtkm/worklet/OrientPointNormals.h

@@ -0,0 +1,369 @@
+//============================================================================
+//  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 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2019 UT-Battelle, LLC.
+//  Copyright 2019 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 vtkm_m_worklet_OrientPointNormals_h
+#define vtkm_m_worklet_OrientPointNormals_h
+
+#include <vtkm/Range.h>
+#include <vtkm/Types.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/ArrayHandleBitField.h>
+#include <vtkm/cont/ArrayHandleConstant.h>
+#include <vtkm/cont/ArrayHandleTransform.h>
+#include <vtkm/cont/ArrayRangeCompute.h>
+#include <vtkm/cont/BitField.h>
+#include <vtkm/cont/Logging.h>
+
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/MaskIndices.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+
+///
+/// Orients normals to point outside of the dataset. This requires a closed
+/// manifold surface or else the behavior is undefined. This requires an
+/// unstructured cellset as input.
+///
+class OrientPointNormals
+{
+  static constexpr vtkm::Id INVALID_ID = -1;
+
+  // Returns true if v1 and v2 are pointing in the same hemisphere.
+  template <typename T>
+  VTKM_EXEC static bool SameDirection(const vtkm::Vec<T, 3>& v1, const vtkm::Vec<T, 3>& v2)
+  {
+    return vtkm::Dot(v1, v2) >= 0;
+  }
+
+  // Ensure that the normal is pointing in the same hemisphere as ref.
+  // Returns true if normal is modified.
+  template <typename T>
+  VTKM_EXEC static bool Align(vtkm::Vec<T, 3>& normal, const vtkm::Vec<T, 3>& ref)
+  {
+    if (!SameDirection(normal, ref))
+    {
+      normal = -normal;
+      return true;
+    }
+    return false;
+  }
+
+public:
+  // Locates starting points for BFS traversal of dataset by finding points
+  // on the dataset boundaries. The normals for these points are corrected by
+  // making them point outside of the dataset, and they are marked as both
+  // active and visited.
+  class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldIn coords,
+                                  FieldInOut normals,
+                                  WholeArrayIn ranges,
+                                  FieldOut activePoints,
+                                  FieldOut visitedPoints,
+                                  FieldOut refPoints);
+    using ExecutionSignature =
+      _6(InputIndex pointId, _1 coord, _2 normal, _3 ranges, _4 activePoints, _5 visitedPoints);
+
+    template <typename CoordT, typename NormalT, typename RangePortal>
+    VTKM_EXEC vtkm::Id operator()(const vtkm::Id pointId,
+                                  const vtkm::Vec<CoordT, 3>& point,
+                                  vtkm::Vec<NormalT, 3>& normal,
+                                  const RangePortal& ranges,
+                                  bool& isActive,
+                                  bool& isVisited) const
+    {
+      for (vtkm::IdComponent dim = 0; dim < 3; ++dim)
+      {
+        const auto& range = ranges.Get(dim);
+        const auto val = static_cast<decltype(range.Min)>(point[dim]);
+        if (val <= range.Min)
+        {
+          vtkm::Vec<NormalT, 3> ref{ static_cast<NormalT>(0) };
+          ref[dim] = static_cast<NormalT>(-1);
+          Align(normal, ref);
+          isActive = true;
+          isVisited = true;
+          return pointId;
+        }
+        else if (val >= range.Max)
+        {
+          vtkm::Vec<NormalT, 3> ref{ static_cast<NormalT>(0) };
+          ref[dim] = static_cast<NormalT>(1);
+          Align(normal, ref);
+          isActive = true;
+          isVisited = true;
+          return pointId;
+        }
+      }
+
+      isActive = false;
+      isVisited = false;
+      return INVALID_ID;
+    }
+  };
+
+  // Traverses the active points (via mask) and marks the connected cells as
+  // active. Set the reference point for all adjacent cells to the current
+  // point.
+  class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint
+  {
+  public:
+    using ControlSignature = void(CellSetIn cellSet,
+                                  // InOut to preserve data on masked indices
+                                  BitFieldInOut activeCells,
+                                  BitFieldInOut visitedCells,
+                                  FieldInOutPoint activePoints);
+    using ExecutionSignature = _4(CellIndices cells, _2 activeCells, _3 visitedCells);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    // Mark all unvisited cells as active:
+    template <typename CellListT, typename ActiveCellsT, typename VisitedCellsT>
+    VTKM_EXEC bool operator()(const CellListT& cells,
+                              ActiveCellsT& activeCells,
+                              VisitedCellsT& visitedCells) const
+    {
+      for (vtkm::IdComponent c = 0; c < cells.GetNumberOfComponents(); ++c)
+      {
+        const vtkm::Id cellId = cells[c];
+        const bool alreadyVisited = visitedCells.CompareAndSwapBitAtomic(cellId, true, false);
+
+        if (!alreadyVisited)
+        { // This thread is first to visit cell
+          activeCells.SetBitAtomic(cellId, true);
+        }
+      }
+
+      // Mark the current point as inactive:
+      return false;
+    }
+  };
+
+  // Traverses the active cells and mark the connected points as active,
+  // propogating the reference pointId.
+  class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell
+  {
+  public:
+    using ControlSignature = void(CellSetIn cellSet,
+                                  BitFieldInOut activePoints,
+                                  BitFieldIn visitedPoints,
+                                  WholeArrayInOut refPoints,
+                                  FieldInOutCell activeCells);
+    using ExecutionSignature = _5(PointIndices points,
+                                  _2 activePoints,
+                                  _3 visitedPoints,
+                                  _4 refPoints);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename PointListT,
+              typename ActivePointsT,
+              typename VisitedPointsT,
+              typename RefPointsT>
+    VTKM_EXEC bool operator()(const PointListT& points,
+                              ActivePointsT& activePoints,
+                              VisitedPointsT& visitedPoints,
+                              RefPointsT& refPoints) const
+    {
+      // Find any point in the cell that has already been visited, and take
+      // its id as the reference for this cell.
+      vtkm::Id refPtId = INVALID_ID;
+      for (vtkm::IdComponent p = 0; p < points.GetNumberOfComponents(); ++p)
+      {
+        const vtkm::Id pointId = points[p];
+        const bool alreadyVisited = visitedPoints.GetBit(pointId);
+        if (alreadyVisited)
+        {
+          refPtId = pointId;
+          break;
+        }
+      }
+
+      // There must be one valid point in each cell:
+      VTKM_ASSERT("Reference point not found." && refPtId != INVALID_ID);
+
+      // Propogate the reference point to other cell members
+      for (vtkm::IdComponent p = 0; p < points.GetNumberOfComponents(); ++p)
+      {
+        const vtkm::Id pointId = points[p];
+
+        // Mark this point as active
+        const bool alreadyVisited = visitedPoints.GetBit(pointId);
+        if (!alreadyVisited)
+        {
+          const bool alreadyActive = activePoints.CompareAndSwapBitAtomic(pointId, true, false);
+          if (!alreadyActive)
+          { // If we're the first thread to mark point active, set ref point:
+            refPoints.Set(pointId, refPtId);
+          }
+        }
+      }
+
+      // Mark current cell as inactive:
+      return false;
+    }
+  };
+
+  // For each point with a refPtId set, ensure that the associated normal is
+  // in the same hemisphere as the reference normal.
+  // This must be done in a separate step from MarkActivePoints since modifying
+  // visitedPoints in that worklet would create race conditions.
+  class WorkletProcessNormals : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldIn refIds,
+                                  WholeArrayInOut normals,
+                                  // InOut to preserve data on masked indices
+                                  BitFieldInOut visitedPoints);
+    using ExecutionSignature = void(InputIndex ptId, _1 refPtId, _2 normals, _3 visitedPoints);
+
+    using MaskType = vtkm::worklet::MaskIndices;
+
+    template <typename NormalsPortal, typename VisitedPointsT>
+    VTKM_EXEC void operator()(const vtkm::Id ptId,
+                              const vtkm::Id refPtId,
+                              NormalsPortal& normals,
+                              VisitedPointsT& visitedPoints) const
+    {
+      visitedPoints.SetBitAtomic(ptId, true);
+
+      using Normal = typename NormalsPortal::ValueType;
+      Normal normal = normals.Get(ptId);
+      const Normal ref = normals.Get(refPtId);
+      if (Align(normal, ref))
+      {
+        normals.Set(ptId, normal);
+      }
+    }
+  };
+
+  template <typename CellSetType,
+            typename CoordsCompType,
+            typename CoordsStorageType,
+            typename PointNormalCompType,
+            typename PointNormalStorageType>
+  VTKM_CONT static void Run(
+    const CellSetType& cells,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsCompType, 3>, CoordsStorageType>& coords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<PointNormalCompType, 3>, PointNormalStorageType>&
+      pointNormals)
+  {
+    using RangeType = vtkm::cont::ArrayHandle<vtkm::Range>;
+
+    using MarkSourcePoints = vtkm::worklet::DispatcherMapField<WorkletMarkSourcePoints>;
+    using MarkActiveCells = vtkm::worklet::DispatcherMapTopology<WorkletMarkActiveCells>;
+    using MarkActivePoints = vtkm::worklet::DispatcherMapTopology<WorkletMarkActivePoints>;
+    using ProcessNormals = vtkm::worklet::DispatcherMapField<WorkletProcessNormals>;
+
+    const vtkm::Id numCells = cells.GetNumberOfCells();
+
+    VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf,
+                   "OrientPointNormals worklet (%lld points, %lld cells)",
+                   static_cast<vtkm::Int64>(coords.GetNumberOfValues()),
+                   static_cast<vtkm::Int64>(numCells));
+
+    // active = cells / point to be used in the next worklet invocation mask.
+    vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints
+    auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits);
+
+    vtkm::cont::BitField activeCellBits;
+    vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells);
+    auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits);
+
+    // visited = cells / points that have been corrected.
+    vtkm::cont::BitField visitedPointBits; // Initialized by MarkSourcePoints
+    auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits);
+
+    vtkm::cont::BitField visitedCellBits;
+    vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells);
+    auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits);
+
+    vtkm::cont::ArrayHandle<vtkm::Id> mask; // Allocated as needed
+
+    // For each point, store a reference alignment point. Allocated by
+    // MarkSourcePoints.
+    vtkm::cont::ArrayHandle<vtkm::Id> refPoints;
+
+    // 1) Compute range of coords.
+    const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords);
+
+    // 2) Label source points for traversal (use those on a boundary).
+    //    Correct the normals for these points by making them point towards the
+    //    boundary.
+    {
+      MarkSourcePoints dispatcher;
+      dispatcher.Invoke(coords, pointNormals, ranges, activePoints, visitedPoints, refPoints);
+    }
+
+    for (size_t iter = 1;; ++iter)
+    {
+      // 3) Mark unvisited cells adjacent to active points
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActiveCells from " << numActive << " active points.");
+        MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, activeCellBits, visitedCellBits, activePoints);
+      }
+
+      // 4) Mark unvisited points in active cells, using ref point from cell.
+      {
+        vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask);
+        (void)numActive;
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                   "MarkActivePoints from " << numActive << " active cells.");
+        MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(cells, activePointBits, visitedPointBits, refPoints, activeCells);
+      }
+
+      vtkm::Id numActivePoints =
+        vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask);
+
+      if (numActivePoints == 0)
+      { // Done!
+        VTKM_LOG_S(vtkm::cont::LogLevel::Perf, "Iteration " << iter << ": Traversal complete.");
+        break;
+      }
+
+      VTKM_LOG_S(vtkm::cont::LogLevel::Perf,
+                 "Iteration " << iter << ": Processing " << numActivePoints << " normals.");
+
+      // 5) Correct normals for active points.
+      {
+        ProcessNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } };
+        dispatcher.Invoke(refPoints, pointNormals, visitedPointBits);
+      }
+    }
+  }
+};
+}
+} // end namespace vtkm::worklet
+
+
+#endif // vtkm_m_worklet_OrientPointNormals_h

vtkm/worklet/testing/CMakeLists.txt

@@ -45,6 +45,7 @@ set(unit_tests
   UnitTestNDimsEntropy.cxx
   UnitTestNDimsHistogram.cxx
   UnitTestNDimsHistMarginalization.cxx
+  UnitTestOrientNormals.cxx
   UnitTestParticleAdvection.cxx
   UnitTestPointElevation.cxx
   UnitTestPointGradient.cxx

vtkm/worklet/testing/UnitTestOrientNormals.cxx

@@ -0,0 +1,408 @@
+//=============================================================================
+//
+//  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 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2019 UT-Battelle, LLC.
+//  Copyright 2019 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/worklet/OrientNormals.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayCopy.h>
+#include <vtkm/cont/ArrayHandleBitField.h>
+#include <vtkm/cont/ArrayHandleConstant.h>
+#include <vtkm/cont/ArrayHandleVirtual.h>
+#include <vtkm/cont/BitField.h>
+#include <vtkm/cont/CellSet.h>
+#include <vtkm/cont/CellSetSingleType.h>
+#include <vtkm/cont/CoordinateSystem.h>
+#include <vtkm/cont/DataSet.h>
+
+#include <vtkm/filter/MarchingCubes.h>
+#include <vtkm/filter/PolicyBase.h>
+#include <vtkm/filter/SurfaceNormals.h>
+
+#include <vtkm/worklet/WaveletGenerator.h>
+
+#include <vtkm/cont/serial/DeviceAdapterSerial.h>
+
+#include <vtkm/cont/testing/Testing.h>
+
+
+namespace
+{
+
+struct TestPolicy : public vtkm::filter::PolicyBase<TestPolicy>
+{
+  using StructuredCellSetList = vtkm::ListTagBase<vtkm::cont::CellSetStructured<3>>;
+  using UnstructuredCellSetList = vtkm::ListTagBase<vtkm::cont::CellSetSingleType<>>;
+  using AllCellSetList = vtkm::ListTagJoin<StructuredCellSetList, UnstructuredCellSetList>;
+  using FieldTypeList = vtkm::ListTagBase<vtkm::FloatDefault, vtkm::Vec<vtkm::FloatDefault, 3>>;
+};
+
+VTKM_CONT
+vtkm::cont::DataSet CreateDataSet(bool pointNormals, bool cellNormals)
+{
+  vtkm::worklet::WaveletGenerator wavelet;
+  wavelet.SetMinimumExtent({ -25 });
+  wavelet.SetMaximumExtent({ 25 });
+  wavelet.SetFrequency({ 20, 15, 25 });
+  wavelet.SetMagnitude({ 5 });
+  auto dataSet = wavelet.GenerateDataSet();
+
+  // Cut a contour
+  vtkm::filter::MarchingCubes contour;
+  contour.SetActiveField("scalars", vtkm::cont::Field::Association::POINTS);
+  contour.SetNumberOfIsoValues(1);
+  contour.SetIsoValue(192);
+  contour.SetMergeDuplicatePoints(true);
+  contour.SetGenerateNormals(false);
+  dataSet = contour.Execute(dataSet, TestPolicy{});
+
+  vtkm::filter::SurfaceNormals normals;
+  normals.SetGeneratePointNormals(pointNormals);
+  normals.SetGenerateCellNormals(cellNormals);
+  normals.SetPointNormalsName("normals");
+  normals.SetCellNormalsName("normals");
+  normals.SetAutoOrientNormals(false);
+  dataSet = normals.Execute(dataSet, TestPolicy{});
+
+  return dataSet;
+}
+
+struct ValidateNormals
+{
+  using CellSetType = vtkm::cont::CellSetSingleType<>;
+  using NormalType = vtkm::Vec<vtkm::FloatDefault, 3>;
+  using NormalsArrayType = vtkm::cont::ArrayHandleVirtual<NormalType>;
+  using NormalsPortalType = decltype(std::declval<NormalsArrayType>().GetPortalConstControl());
+  using ConnType =
+    decltype(std::declval<CellSetType>().PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{},
+                                                         vtkm::TopologyElementTagPoint{},
+                                                         vtkm::TopologyElementTagCell{}));
+  using RConnType =
+    decltype(std::declval<CellSetType>().PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{},
+                                                         vtkm::TopologyElementTagCell{},
+                                                         vtkm::TopologyElementTagPoint{}));
+  using PointsType =
+    decltype(std::declval<vtkm::cont::CoordinateSystem>().GetData().GetPortalConstControl());
+
+  vtkm::cont::CoordinateSystem Coords;
+  CellSetType Cells;
+
+  PointsType Points;
+  ConnType Conn;
+  RConnType RConn;
+
+  NormalsArrayType PointNormalsArray;
+  NormalsPortalType PointNormals;
+  NormalsArrayType CellNormalsArray;
+  NormalsPortalType CellNormals;
+
+  vtkm::cont::BitField VisitedCellsField;
+  vtkm::cont::BitField VisitedPointsField;
+  vtkm::cont::BitField::PortalControl VisitedCells;
+  vtkm::cont::BitField::PortalControl VisitedPoints;
+
+  bool CheckPoints;
+  bool CheckCells;
+
+  VTKM_CONT
+  static void Run(vtkm::cont::DataSet& dataset,
+                  bool checkPoints,
+                  bool checkCells,
+                  const std::string& normalsName = "normals")
+  {
+    vtkm::cont::Field pointNormals;
+    vtkm::cont::Field cellNormals;
+
+    if (checkPoints)
+    {
+      pointNormals = dataset.GetField(normalsName, vtkm::cont::Field::Association::POINTS);
+    }
+    if (checkCells)
+    {
+      cellNormals = dataset.GetField(normalsName, vtkm::cont::Field::Association::CELL_SET);
+    }
+
+    ValidateNormals obj{ dataset, checkPoints, checkCells, pointNormals, cellNormals };
+    obj.Validate();
+  }
+
+  VTKM_CONT
+  ValidateNormals(const vtkm::cont::DataSet& dataset,
+                  bool checkPoints,
+                  bool checkCells,
+                  const vtkm::cont::Field& pointNormalsField,
+                  const vtkm::cont::Field& cellNormalsField)
+    : Coords{ dataset.GetCoordinateSystem() }
+    , Cells{ dataset.GetCellSet().Cast<CellSetType>() }
+    , Points{ this->Coords.GetData().GetPortalConstControl() }
+    , CheckPoints(checkPoints)
+    , CheckCells(checkCells)
+  {
+    // FIXME This would be much simplier if we had a GetPointCells() method on
+    // cell sets....
+    // Build the connectivity table on any device, then get a portal for serial
+    // so we can do lookups on the CPU.
+    this->Cells.GetConnectivityArray(vtkm::TopologyElementTagPoint{},
+                                     vtkm::TopologyElementTagCell{});
+    this->Cells.GetConnectivityArray(vtkm::TopologyElementTagCell{},
+                                     vtkm::TopologyElementTagPoint{});
+    this->Conn = this->Cells.PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{},
+                                             vtkm::TopologyElementTagPoint{},
+                                             vtkm::TopologyElementTagCell{});
+    this->RConn = this->Cells.PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{},
+                                              vtkm::TopologyElementTagCell{},
+                                              vtkm::TopologyElementTagPoint{});
+
+    if (this->CheckPoints)
+    {
+      this->PointNormalsArray = pointNormalsField.GetData().AsVirtual<NormalType>();
+      this->PointNormals = this->PointNormalsArray.GetPortalConstControl();
+    }
+    if (this->CheckCells)
+    {
+      this->CellNormalsArray = cellNormalsField.GetData().AsVirtual<NormalType>();
+      this->CellNormals = this->CellNormalsArray.GetPortalConstControl();
+    }
+  }
+
+  VTKM_CONT
+  void Validate()
+  {
+    // Locate a point with the minimum x coordinate:
+    const vtkm::Id startPoint = [&]() -> vtkm::Id {
+      const vtkm::Float64 xMin = this->Coords.GetBounds().X.Min;
+      const auto points = this->Coords.GetData().GetPortalConstControl();
+      const vtkm::Id numPoints = points.GetNumberOfValues();
+      vtkm::Id resultIdx = -1;
+      for (vtkm::Id pointIdx = 0; pointIdx < numPoints; ++pointIdx)
+      {
+        const auto point = points.Get(pointIdx);
+        if (static_cast<double>(point[0]) <= xMin)
+        {
+          resultIdx = pointIdx;
+          break;
+        }
+      }
+      if (resultIdx < 0)
+      {
+        throw vtkm::cont::ErrorBadValue("Minimum point not found!");
+      }
+
+      return resultIdx;
+    }();
+
+    // Start recursive validation.
+    this->Prepare();
+    this->ValidateImpl(startPoint, NormalType{ -1, 0, 0 });
+
+    vtkm::Id numPoints = this->Points.GetNumberOfValues();
+    vtkm::Id numCells = this->Cells.GetNumberOfCells();
+    vtkm::Id numVisitedPoints = vtkm::cont::Algorithm::CountSetBits(this->VisitedPointsField);
+    vtkm::Id numVisitedCells = vtkm::cont::Algorithm::CountSetBits(this->VisitedCellsField);
+    if (numPoints != numVisitedPoints)
+    {
+      throw vtkm::cont::ErrorBadValue("Unvisited point!");
+    }
+    if (numCells != numVisitedCells)
+    {
+      throw vtkm::cont::ErrorBadValue("Unvisited cell!");
+    }
+  }
+
+private:
+  static bool SameHemisphere(const NormalType& a, const NormalType& b)
+  {
+    return vtkm::Dot(a, b) >= 0;
+  }
+
+  void Prepare()
+  {
+    vtkm::cont::Algorithm::Fill(this->VisitedPointsField, false, this->Coords.GetNumberOfPoints());
+    vtkm::cont::Algorithm::Fill(this->VisitedCellsField, false, this->Cells.GetNumberOfCells());
+
+    this->VisitedPoints = this->VisitedPointsField.GetPortalControl();
+    this->VisitedCells = this->VisitedCellsField.GetPortalControl();
+  }
+
+  void ValidateImpl(vtkm::Id startPtIdx, const NormalType& startRefNormal)
+  {
+    using Entry = vtkm::Pair<vtkm::Id, NormalType>;
+    std::vector<Entry> queue;
+    queue.emplace_back(startPtIdx, startRefNormal);
+    this->VisitedPoints.SetBit(startPtIdx, true);
+
+    while (!queue.empty())
+    {
+      const vtkm::Id curPtIdx = queue.back().first;
+      NormalType refNormal = queue.back().second;
+      queue.pop_back();
+
+      if (this->CheckPoints)
+      {
+        const NormalType curNormal = this->PointNormals.Get(curPtIdx);
+        if (!this->SameHemisphere(curNormal, refNormal))
+        {
+          const auto coord = this->Points.Get(curPtIdx);
+          std::cerr << "BadPointNormal PtId: " << curPtIdx << "\n\t"
+                    << "- Normal: {" << curNormal[0] << ", " << curNormal[1] << ", " << curNormal[2]
+                    << "}\n\t"
+                    << "- Reference: {" << refNormal[0] << ", " << refNormal[1] << ", "
+                    << refNormal[2] << "}\n\t"
+                    << "- Coord: {" << coord[0] << ", " << coord[1] << ", " << coord[2] << "}\n";
+          throw vtkm::cont::ErrorBadValue("Bad point normal found!");
+        }
+        refNormal = curNormal;
+      }
+
+      // Lookup and visit neighbor cells:
+      const auto neighborCells = this->RConn.GetIndices(curPtIdx);
+      const auto numNeighborCells = neighborCells.GetNumberOfComponents();
+      for (vtkm::IdComponent nCellIdx = 0; nCellIdx < numNeighborCells; ++nCellIdx)
+      {
+        const vtkm::Id curCellIdx = neighborCells[nCellIdx];
+
+        // Skip this cell if already visited:
+        if (this->VisitedCells.GetBit(curCellIdx))
+        {
+          continue;
+        }
+        this->VisitedCells.SetBit(curCellIdx, true);
+
+        if (this->CheckCells)
+        {
+          const NormalType curNormal = this->CellNormals.Get(curCellIdx);
+          if (!this->SameHemisphere(curNormal, refNormal))
+          {
+            std::cerr << "BadCellNormal: CellId: " << curCellIdx << "\n\t"
+                      << "- Normal: {" << curNormal[0] << ", " << curNormal[1] << ", "
+                      << curNormal[2] << "}\n\t"
+                      << "- Reference: {" << refNormal[0] << ", " << refNormal[1] << ", "
+                      << refNormal[2] << "}\n";
+            throw vtkm::cont::ErrorBadValue("Bad cell normal found!");
+          }
+          refNormal = curNormal;
+        }
+
+        // Lookup and visit points in this cell:
+        const auto neighborPoints = this->Conn.GetIndices(curCellIdx);
+        const auto numNeighborPoints = neighborPoints.GetNumberOfComponents();
+        for (vtkm::IdComponent nPtIdx = 0; nPtIdx < numNeighborPoints; ++nPtIdx)
+        {
+          const vtkm::Id nextPtIdx = neighborPoints[nPtIdx];
+
+          // Skip if already visited:
+          if (this->VisitedPoints.GetBit(nextPtIdx))
+          {
+            continue;
+          }
+
+          // Otherwise, queue next point using current normal as reference:
+          queue.emplace_back(nextPtIdx, refNormal);
+          this->VisitedPoints.SetBit(nextPtIdx, true);
+        }
+      }
+    }
+  }
+};
+
+VTKM_CONT
+void TestOrientNormals(bool testPoints, bool testCells)
+{
+  using NormalArrayT = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
+
+  auto dataset = CreateDataSet(testPoints, testCells);
+
+  // Check that the input actually has bad normals:
+  const bool inputValid = [&]() -> bool {
+    try
+    {
+      std::cerr << "Expecting to throw a BadNormal exception...\n";
+      ValidateNormals::Run(dataset, testPoints, testCells);
+      return true; // Dataset is already oriented
+    }
+    catch (vtkm::cont::ErrorBadValue&)
+    {
+      std::cerr << "Expected exception caught! All is well. "
+                   "Future exceptions are errors.\n";
+      return false; // Dataset is unoriented
+    }
+  }();
+
+  if (inputValid)
+  {
+    throw vtkm::cont::ErrorBadValue("Error: Input doesn't have bad normals.");
+  }
+
+  // modify normals in place
+  const auto coords = dataset.GetCoordinateSystem().GetData();
+  const auto cells = dataset.GetCellSet();
+  if (testPoints && testCells)
+  {
+    std::cerr << "Testing point and cell normals...\n";
+
+    const auto pointNormalField =
+      dataset.GetField("normals", vtkm::cont::Field::Association::POINTS);
+    const auto cellNormalField =
+      dataset.GetField("normals", vtkm::cont::Field::Association::CELL_SET);
+    auto pointNormals = pointNormalField.GetData().Cast<NormalArrayT>();
+    auto cellNormals = cellNormalField.GetData().Cast<NormalArrayT>();
+
+    vtkm::worklet::OrientNormals::RunPointAndCellNormals(cells, coords, pointNormals, cellNormals);
+  }
+  else if (testPoints)
+  {
+    std::cerr << "Testing point normals...\n";
+    const auto pointNormalField =
+      dataset.GetField("normals", vtkm::cont::Field::Association::POINTS);
+    auto pointNormals = pointNormalField.GetData().Cast<NormalArrayT>();
+
+    vtkm::worklet::OrientNormals::RunPointNormals(cells, coords, pointNormals);
+  }
+  else if (testCells)
+  {
+    std::cerr << "Testing cell normals...\n";
+    const auto cellNormalField =
+      dataset.GetField("normals", vtkm::cont::Field::Association::CELL_SET);
+    auto cellNormals = cellNormalField.GetData().Cast<NormalArrayT>();
+
+    vtkm::worklet::OrientNormals::RunCellNormals(cells, coords, cellNormals);
+  }
+  else
+  {
+    throw "Nothing tested...";
+  }
+
+  ValidateNormals::Run(dataset, testPoints, testCells);
+}
+
+void DoTest()
+{
+  TestOrientNormals(true, false);
+  TestOrientNormals(false, true);
+  TestOrientNormals(true, true);
+}
+
+} // end anon namespace
+
+int UnitTestOrientNormals(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(DoTest, argc, argv);
+}