diff --git a/vtkm/rendering/CMakeLists.txt b/vtkm/rendering/CMakeLists.txt
index 03ee5a037..03dca45f9 100644
--- a/vtkm/rendering/CMakeLists.txt
+++ b/vtkm/rendering/CMakeLists.txt
@@ -178,6 +178,7 @@ set(device_sources
   raytracing/SphereExtractor.cxx
   raytracing/SphereIntersector.cxx
   raytracing/TriangleExtractor.cxx
+  raytracing/TriangleIntersector.cxx
   raytracing/VolumeRendererStructured.cxx
   )
 
diff --git a/vtkm/rendering/ConnectivityProxy.cxx b/vtkm/rendering/ConnectivityProxy.cxx
index 245af777a..661abb9db 100644
--- a/vtkm/rendering/ConnectivityProxy.cxx
+++ b/vtkm/rendering/ConnectivityProxy.cxx
@@ -68,6 +68,9 @@ public:
 
   ~InternalsType() {}
 
+  VTKM_CONT
+  void SetUnitScalar(vtkm::Float32 unitScalar) { Tracer.SetUnitScalar(unitScalar); }
+
   void SetSampleDistance(const vtkm::Float32& distance)
   {
     if (Mode != VOLUME_MODE)
@@ -119,7 +122,12 @@ public:
   vtkm::Bounds GetSpatialBounds() const { return SpatialBounds; }
 
   VTKM_CONT
-  vtkm::Range GetScalarRange() const { return ScalarRange; }
+  vtkm::Range GetScalarFieldRange()
+  {
+    const vtkm::cont::ArrayHandle<vtkm::Range> range = this->ScalarField.GetRange();
+    ScalarRange = range.GetPortalConstControl().Get(0);
+    return ScalarRange;
+  }
 
   VTKM_CONT
   void SetScalarRange(const vtkm::Range& range) { ScalarRange = range; }
@@ -141,7 +149,7 @@ public:
                            this->EmissionField);
     }
 
-    Tracer.Trace(rays);
+    Tracer.FullTrace(rays);
   }
 
   VTKM_CONT
@@ -160,7 +168,46 @@ public:
                            this->EmissionField);
     }
 
-    Tracer.Trace(rays);
+    Tracer.FullTrace(rays);
+  }
+
+  VTKM_CONT
+  PartialVector64 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays)
+  {
+
+    if (Mode == VOLUME_MODE)
+    {
+      Tracer.SetVolumeData(this->ScalarField, this->ScalarRange, this->Cells, this->Coords);
+    }
+    else
+    {
+      Tracer.SetEnergyData(this->ScalarField,
+                           rays.Buffers.at(0).GetNumChannels(),
+                           this->Cells,
+                           this->Coords,
+                           this->EmissionField);
+    }
+
+    return Tracer.PartialTrace(rays);
+  }
+
+  VTKM_CONT
+  PartialVector32 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
+  {
+    if (Mode == VOLUME_MODE)
+    {
+      Tracer.SetVolumeData(this->ScalarField, this->ScalarRange, this->Cells, this->Coords);
+    }
+    else
+    {
+      Tracer.SetEnergyData(this->ScalarField,
+                           rays.Buffers.at(0).GetNumChannels(),
+                           this->Cells,
+                           this->Coords,
+                           this->EmissionField);
+    }
+
+    return Tracer.PartialTrace(rays);
   }
 
   VTKM_CONT
@@ -188,7 +235,7 @@ public:
       throw vtkm::cont::ErrorBadValue("ENERGY MODE Not implemented for this use case\n");
     }
 
-    Tracer.Trace(rays);
+    Tracer.FullTrace(rays);
 
     canvas->WriteToCanvas(rays, rays.Buffers.at(0).Buffer, camera);
     if (CompositeBackground)
@@ -266,9 +313,9 @@ vtkm::Bounds ConnectivityProxy::GetSpatialBounds()
 }
 
 VTKM_CONT
-vtkm::Range ConnectivityProxy::GetScalarRange()
+vtkm::Range ConnectivityProxy::GetScalarFieldRange()
 {
-  return Internals->GetScalarRange();
+  return Internals->GetScalarFieldRange();
 }
 
 VTKM_CONT
@@ -301,6 +348,27 @@ void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& r
   logger->CloseLogEntry(-1.0);
 }
 
+VTKM_CONT
+PartialVector32 ConnectivityProxy::PartialTrace(
+  vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
+{
+  raytracing::Logger* logger = raytracing::Logger::GetInstance();
+  logger->OpenLogEntry("connectivity_trace_32");
+  if (Internals->GetRenderMode() == VOLUME_MODE)
+  {
+    logger->AddLogData("volume_mode", "true");
+  }
+  else
+  {
+    logger->AddLogData("volume_mode", "false");
+  }
+
+  PartialVector32 res = Internals->PartialTrace(rays);
+
+  logger->CloseLogEntry(-1.0);
+  return res;
+}
+
 VTKM_CONT
 void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
 {
@@ -320,6 +388,27 @@ void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& r
   logger->CloseLogEntry(-1.0);
 }
 
+VTKM_CONT
+PartialVector64 ConnectivityProxy::PartialTrace(
+  vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays)
+{
+  raytracing::Logger* logger = raytracing::Logger::GetInstance();
+  logger->OpenLogEntry("connectivity_trace_64");
+  if (Internals->GetRenderMode() == VOLUME_MODE)
+  {
+    logger->AddLogData("volume_mode", "true");
+  }
+  else
+  {
+    logger->AddLogData("volume_mode", "false");
+  }
+
+  PartialVector64 res = Internals->PartialTrace(rays);
+
+  logger->CloseLogEntry(-1.0);
+  return res;
+}
+
 VTKM_CONT
 void ConnectivityProxy::Trace(const vtkm::rendering::Camera& camera,
                               vtkm::rendering::CanvasRayTracer* canvas)
@@ -332,10 +421,17 @@ void ConnectivityProxy::Trace(const vtkm::rendering::Camera& camera,
 
   logger->CloseLogEntry(-1.0);
 }
+
 VTKM_CONT
 void ConnectivityProxy::SetDebugPrints(bool on)
 {
   Internals->SetDebugPrints(on);
 }
+
+VTKM_CONT
+void ConnectivityProxy::SetUnitScalar(vtkm::Float32 unitScalar)
+{
+  Internals->SetUnitScalar(unitScalar);
+}
 }
 } // namespace vtkm::rendering
diff --git a/vtkm/rendering/ConnectivityProxy.h b/vtkm/rendering/ConnectivityProxy.h
index fec9634a7..ba73d9c94 100644
--- a/vtkm/rendering/ConnectivityProxy.h
+++ b/vtkm/rendering/ConnectivityProxy.h
@@ -26,6 +26,7 @@
 #include <vtkm/rendering/Mapper.h>
 #include <vtkm/rendering/View.h>
 #include <vtkm/rendering/raytracing/Camera.h>
+#include <vtkm/rendering/raytracing/PartialComposite.h>
 #include <vtkm/rendering/raytracing/Ray.h>
 
 namespace vtkm
@@ -33,6 +34,9 @@ namespace vtkm
 namespace rendering
 {
 
+using PartialVector64 = std::vector<vtkm::rendering::raytracing::PartialComposite<vtkm::Float64>>;
+using PartialVector32 = std::vector<vtkm::rendering::raytracing::PartialComposite<vtkm::Float32>>;
+
 class VTKM_RENDERING_EXPORT ConnectivityProxy
 {
 public:
@@ -57,14 +61,18 @@ public:
   void SetColorMap(vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>& colormap);
   void SetCompositeBackground(bool on);
   void SetDebugPrints(bool on);
+  void SetUnitScalar(vtkm::Float32 unitScalar);
 
   vtkm::Bounds GetSpatialBounds();
-  vtkm::Range GetScalarRange();
+  vtkm::Range GetScalarFieldRange();
 
   void Trace(const vtkm::rendering::Camera& camera, vtkm::rendering::CanvasRayTracer* canvas);
   void Trace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays);
   void Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays);
 
+  PartialVector64 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays);
+  PartialVector32 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays);
+
 protected:
   struct InternalsType;
   struct BoundsFunctor;
diff --git a/vtkm/rendering/MapperCylinder.cxx b/vtkm/rendering/MapperCylinder.cxx
index 20ab009f2..0a176a314 100644
--- a/vtkm/rendering/MapperCylinder.cxx
+++ b/vtkm/rendering/MapperCylinder.cxx
@@ -29,12 +29,8 @@
 #include <vtkm/rendering/raytracing/CylinderExtractor.h>
 #include <vtkm/rendering/raytracing/CylinderIntersector.h>
 #include <vtkm/rendering/raytracing/Logger.h>
-#include <vtkm/rendering/raytracing/QuadExtractor.h>
-#include <vtkm/rendering/raytracing/QuadIntersector.h>
 #include <vtkm/rendering/raytracing/RayOperations.h>
 #include <vtkm/rendering/raytracing/RayTracer.h>
-#include <vtkm/rendering/raytracing/SphereExtractor.h>
-#include <vtkm/rendering/raytracing/SphereIntersector.h>
 #include <vtkm/rendering/raytracing/Worklets.h>
 
 namespace vtkm
diff --git a/vtkm/rendering/MapperPoint.cxx b/vtkm/rendering/MapperPoint.cxx
index 169120ef6..fa401d285 100644
--- a/vtkm/rendering/MapperPoint.cxx
+++ b/vtkm/rendering/MapperPoint.cxx
@@ -21,7 +21,6 @@
 #include <vtkm/rendering/MapperPoint.h>
 
 #include <vtkm/cont/Timer.h>
-#include <vtkm/cont/TryExecute.h>
 
 #include <vtkm/rendering/CanvasRayTracer.h>
 #include <vtkm/rendering/internal/RunTriangulator.h>
@@ -31,7 +30,6 @@
 #include <vtkm/rendering/raytracing/RayTracer.h>
 #include <vtkm/rendering/raytracing/SphereExtractor.h>
 #include <vtkm/rendering/raytracing/SphereIntersector.h>
-#include <vtkm/rendering/raytracing/TriangleExtractor.h>
 
 namespace vtkm
 {
diff --git a/vtkm/rendering/MapperQuad.cxx b/vtkm/rendering/MapperQuad.cxx
index 4c3f62192..f96c47b93 100644
--- a/vtkm/rendering/MapperQuad.cxx
+++ b/vtkm/rendering/MapperQuad.cxx
@@ -26,15 +26,11 @@
 #include <vtkm/rendering/CanvasRayTracer.h>
 #include <vtkm/rendering/Cylinderizer.h>
 #include <vtkm/rendering/raytracing/Camera.h>
-#include <vtkm/rendering/raytracing/CylinderExtractor.h>
-#include <vtkm/rendering/raytracing/CylinderIntersector.h>
 #include <vtkm/rendering/raytracing/Logger.h>
 #include <vtkm/rendering/raytracing/QuadExtractor.h>
 #include <vtkm/rendering/raytracing/QuadIntersector.h>
 #include <vtkm/rendering/raytracing/RayOperations.h>
 #include <vtkm/rendering/raytracing/RayTracer.h>
-#include <vtkm/rendering/raytracing/SphereExtractor.h>
-#include <vtkm/rendering/raytracing/SphereIntersector.h>
 
 namespace vtkm
 {
diff --git a/vtkm/rendering/raytracing/BVHTraverser.h b/vtkm/rendering/raytracing/BVHTraverser.h
index e3b38c457..9b1dc6332 100644
--- a/vtkm/rendering/raytracing/BVHTraverser.h
+++ b/vtkm/rendering/raytracing/BVHTraverser.h
@@ -78,28 +78,13 @@ VTKM_EXEC inline bool IntersectAABB(const BVHPortalType& bvh,
   return (min0 > min1);
 }
 
-template <template <typename> class LeafIntersectorType>
 class BVHTraverser
 {
 public:
-  template <typename Device>
   class Intersector : public vtkm::worklet::WorkletMapField
   {
-  public:
-    typedef typename vtkm::cont::ArrayHandle<Vec<vtkm::Float32, 4>> Float4ArrayHandle;
-    typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Int32, 2>> Int2Handle;
-    typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Id4Handle;
-    typedef typename vtkm::cont::ArrayHandle<vtkm::Id> IdHandle;
-    typedef typename Float4ArrayHandle::ExecutionTypes<Device>::PortalConst Float4ArrayPortal;
-    typedef typename Int2Handle::ExecutionTypes<Device>::PortalConst Int2ArrayPortal;
-    typedef typename IdHandle::ExecutionTypes<Device>::PortalConst IdArrayPortal;
-    typedef typename Id4Handle::ExecutionTypes<Device>::PortalConst Id4ArrayPortal;
-
   private:
     bool Occlusion;
-    Float4ArrayPortal FlatBVH;
-    IdArrayPortal Leafs;
-    LeafIntersectorType<Device> LeafIntersector;
 
     VTKM_EXEC
     inline vtkm::Float32 rcp(vtkm::Float32 f) const { return 1.0f / f; }
@@ -118,11 +103,8 @@ public:
 
   public:
     VTKM_CONT
-    Intersector(bool occlusion, LinearBVH& bvh, LeafIntersectorType<Device>& leafIntersector)
+    Intersector(bool occlusion)
       : Occlusion(occlusion)
-      , FlatBVH(bvh.FlatBVH.PrepareForInput(Device()))
-      , Leafs(bvh.Leafs.PrepareForInput(Device()))
-      , LeafIntersector(leafIntersector)
     {
     }
     using ControlSignature = void(FieldIn<>,
@@ -133,11 +115,18 @@ public:
                                   FieldOut<>,
                                   FieldOut<>,
                                   FieldOut<>,
-                                  WholeArrayIn<Vec3RenderingTypes>);
-    using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9);
+                                  WholeArrayIn<Vec3RenderingTypes>,
+                                  ExecObject leafIntersector,
+                                  WholeArrayIn<>,
+                                  WholeArrayIn<>);
+    using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12);
 
 
-    template <typename PointPortalType, typename Precision>
+    template <typename PointPortalType,
+              typename Precision,
+              typename LeafType,
+              typename InnerNodePortalType,
+              typename LeafPortalType>
     VTKM_EXEC void operator()(const vtkm::Vec<Precision, 3>& dir,
                               const vtkm::Vec<Precision, 3>& origin,
                               Precision& distance,
@@ -146,7 +135,10 @@ public:
                               Precision& minU,
                               Precision& minV,
                               vtkm::Id& hitIndex,
-                              const PointPortalType& points) const
+                              const PointPortalType& points,
+                              LeafType& leafIntersector,
+                              const InnerNodePortalType& flatBVH,
+                              const LeafPortalType& leafs) const
     {
       Precision closestDistance = maxDistance;
       distance = maxDistance;
@@ -174,7 +166,7 @@ public:
 
 
           bool hitLeftChild, hitRightChild;
-          bool rightCloser = IntersectAABB(FlatBVH,
+          bool rightCloser = IntersectAABB(flatBVH,
                                            currentNode,
                                            originDir,
                                            invDir,
@@ -191,7 +183,7 @@ public:
           else
           {
             vtkm::Vec<vtkm::Float32, 4> children =
-              FlatBVH.Get(currentNode + 3); //Children.Get(currentNode);
+              flatBVH.Get(currentNode + 3); //Children.Get(currentNode);
             vtkm::Int32 leftChild;
             memcpy(&leftChild, &children[0], 4);
             vtkm::Int32 rightChild;
@@ -217,7 +209,7 @@ public:
         if (currentNode < 0 && currentNode != barrier) //check register usage
         {
           currentNode = -currentNode - 1; //swap the neg address
-          LeafIntersector.IntersectLeaf(currentNode,
+          leafIntersector.IntersectLeaf(currentNode,
                                         origin,
                                         dir,
                                         points,
@@ -225,7 +217,7 @@ public:
                                         closestDistance,
                                         minU,
                                         minV,
-                                        Leafs,
+                                        leafs,
                                         minDistance);
           currentNode = todo[stackptr];
           stackptr--;
@@ -239,16 +231,13 @@ public:
   };
 
 
-  template <typename Precision, typename Device>
+  template <typename Precision, typename LeafIntersectorType>
   VTKM_CONT void IntersectRays(Ray<Precision>& rays,
                                LinearBVH& bvh,
-                               LeafIntersectorType<Device> leafIntersector,
-                               vtkm::cont::CoordinateSystem& coordsHandle,
-                               Device vtkmNotUsed(Device))
+                               LeafIntersectorType& leafIntersector,
+                               vtkm::cont::CoordinateSystem& coordsHandle)
   {
-    vtkm::worklet::DispatcherMapField<Intersector<Device>> intersectDispatch(
-      Intersector<Device>(false, bvh, leafIntersector));
-    intersectDispatch.SetDevice(Device());
+    vtkm::worklet::DispatcherMapField<Intersector> intersectDispatch(Intersector(false));
     intersectDispatch.Invoke(rays.Dir,
                              rays.Origin,
                              rays.Distance,
@@ -257,7 +246,10 @@ public:
                              rays.U,
                              rays.V,
                              rays.HitIdx,
-                             coordsHandle);
+                             coordsHandle,
+                             leafIntersector,
+                             bvh.FlatBVH,
+                             bvh.Leafs);
   }
 }; // BVHTraverser
 #undef END_FLAG
diff --git a/vtkm/rendering/raytracing/CMakeLists.txt b/vtkm/rendering/raytracing/CMakeLists.txt
index 822b2dfca..67da5b672 100644
--- a/vtkm/rendering/raytracing/CMakeLists.txt
+++ b/vtkm/rendering/raytracing/CMakeLists.txt
@@ -35,6 +35,7 @@ set(headers
   MeshConnectivityContainers.h
   MeshConnectivityBase.h
   MortonCodes.h
+  PartialComposite.h
   QuadExtractor.h
   QuadIntersector.h
   Ray.h
@@ -46,6 +47,7 @@ set(headers
   SphereExtractor.h
   SphereIntersector.h
   TriangleExtractor.h
+  TriangleIntersections.h
   TriangleIntersector.h
   VolumeRendererStructured.h
   Worklets.h
diff --git a/vtkm/rendering/raytracing/CellIntersector.h b/vtkm/rendering/raytracing/CellIntersector.h
index ab7ba90fb..908945d68 100644
--- a/vtkm/rendering/raytracing/CellIntersector.h
+++ b/vtkm/rendering/raytracing/CellIntersector.h
@@ -20,8 +20,9 @@
 #ifndef vtk_m_rendering_raytracing_Cell_Intersector_h
 #define vtk_m_rendering_raytracing_Cell_Intersector_h
 
+#include <vtkm/CellShape.h>
 #include <vtkm/rendering/raytracing/CellTables.h>
-#include <vtkm/rendering/raytracing/TriangleIntersector.h>
+#include <vtkm/rendering/raytracing/TriangleIntersections.h>
 
 namespace vtkm
 {
diff --git a/vtkm/rendering/raytracing/ChannelBuffer.cxx b/vtkm/rendering/raytracing/ChannelBuffer.cxx
index b7bc2ba34..4149a388b 100644
--- a/vtkm/rendering/raytracing/ChannelBuffer.cxx
+++ b/vtkm/rendering/raytracing/ChannelBuffer.cxx
@@ -511,6 +511,15 @@ void ChannelBuffer<Precision>::SetNumChannels(const vtkm::Int32 numChannels)
   ResizeChannelFunctor<Precision> functor(this, numChannels);
   vtkm::cont::TryExecute(functor);
 }
+
+template <typename Precision>
+ChannelBuffer<Precision> ChannelBuffer<Precision>::Copy()
+{
+  ChannelBuffer res(NumChannels, Size);
+  vtkm::cont::Algorithm::Copy(this->Buffer, res.Buffer);
+  return res;
+}
+
 // Instantiate supported types
 template class ChannelBuffer<vtkm::Float32>;
 template class ChannelBuffer<vtkm::Float64>;
diff --git a/vtkm/rendering/raytracing/ChannelBuffer.h b/vtkm/rendering/raytracing/ChannelBuffer.h
index e8d13ea3b..510debb8b 100644
--- a/vtkm/rendering/raytracing/ChannelBuffer.h
+++ b/vtkm/rendering/raytracing/ChannelBuffer.h
@@ -83,6 +83,8 @@ public:
                                         const vtkm::Id outputSize,
                                         Precision initValue = 1.f);
 
+  ChannelBuffer<Precision> Copy();
+
   void InitConst(const Precision value);
   void InitChannels(const vtkm::cont::ArrayHandle<Precision>& signature);
   void Normalize(bool invert);
diff --git a/vtkm/rendering/raytracing/ChannelBufferOperations.h b/vtkm/rendering/raytracing/ChannelBufferOperations.h
index 05a193f8c..b2cd5acc4 100644
--- a/vtkm/rendering/raytracing/ChannelBufferOperations.h
+++ b/vtkm/rendering/raytracing/ChannelBufferOperations.h
@@ -103,23 +103,21 @@ public:
 class ChannelBufferOperations
 {
 public:
-  template <typename Device, typename Precision>
+  template <typename Precision>
   static void Compact(ChannelBuffer<Precision>& buffer,
                       vtkm::cont::ArrayHandle<UInt8>& masks,
-                      const vtkm::Id& newSize,
-                      Device)
+                      const vtkm::Id& newSize)
   {
     vtkm::cont::ArrayHandle<vtkm::Id> offsets;
-    offsets.PrepareForOutput(buffer.Size, Device());
+    offsets.Allocate(buffer.Size);
     vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
-    vtkm::cont::DeviceAdapterAlgorithm<Device>::ScanExclusive(castedMasks, offsets);
+    vtkm::cont::Algorithm::ScanExclusive(castedMasks, offsets);
 
     vtkm::cont::ArrayHandle<Precision> compactedBuffer;
-    compactedBuffer.PrepareForOutput(newSize * buffer.NumChannels, Device());
+    compactedBuffer.Allocate(newSize * buffer.NumChannels);
 
     vtkm::worklet::DispatcherMapField<detail::CompactBuffer> dispatcher(
       detail::CompactBuffer(buffer.NumChannels));
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(masks, buffer.Buffer, offsets, compactedBuffer);
     buffer.Buffer = compactedBuffer;
     buffer.Size = newSize;
diff --git a/vtkm/rendering/raytracing/ConnectivityTracer.cxx b/vtkm/rendering/raytracing/ConnectivityTracer.cxx
index 8718fe7c8..c879004df 100644
--- a/vtkm/rendering/raytracing/ConnectivityTracer.cxx
+++ b/vtkm/rendering/raytracing/ConnectivityTracer.cxx
@@ -17,11 +17,39 @@
 //  Laboratory (LANL), the U.S. Government retains certain rights in
 //  this software.
 //============================================================================
-#define vtk_m_rendering_raytracing_ConnectivityTracer_cxx
-
 #include <vtkm/rendering/raytracing/ConnectivityTracer.h>
-#include <vtkm/rendering/raytracing/ConnectivityTracer.hxx>
+
+#include <vtkm/cont/DeviceAdapterAlgorithm.h>
+#include <vtkm/cont/ErrorBadValue.h>
+#include <vtkm/cont/Timer.h>
+#include <vtkm/cont/TryExecute.h>
+#include <vtkm/cont/internal/DeviceAdapterListHelpers.h>
+
+#include <vtkm/rendering/raytracing/Camera.h>
+#include <vtkm/rendering/raytracing/CellIntersector.h>
+#include <vtkm/rendering/raytracing/CellSampler.h>
+#include <vtkm/rendering/raytracing/CellTables.h>
+#include <vtkm/rendering/raytracing/MeshConnectivityBase.h>
 #include <vtkm/rendering/raytracing/MeshConnectivityBuilder.h>
+#include <vtkm/rendering/raytracing/Ray.h>
+#include <vtkm/rendering/raytracing/RayOperations.h>
+#include <vtkm/rendering/raytracing/RayTracingTypeDefs.h>
+#include <vtkm/rendering/raytracing/Worklets.h>
+
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+
+#include <iomanip>
+
+#ifndef CELL_SHAPE_ZOO
+#define CELL_SHAPE_ZOO 255
+#endif
+
+#ifndef CELL_SHAPE_STRUCTURED
+#define CELL_SHAPE_STRUCTURED 254
+#endif
 
 namespace vtkm
 {
@@ -29,33 +57,86 @@ namespace rendering
 {
 namespace raytracing
 {
-
 namespace detail
 {
-struct RenderFunctor
+
+template <typename FloatType>
+void RayTracking<FloatType>::Compact(vtkm::cont::ArrayHandle<FloatType>& compactedDistances,
+                                     vtkm::cont::ArrayHandle<UInt8>& masks)
 {
-  template <typename Device, typename Tracer, typename Rays>
-  bool operator()(Device device, Tracer&& tracer, Rays&& rays) const
+  //
+  // These distances are stored in the rays, and it has
+  // already been compacted.
+  //
+  CurrentDistance = compactedDistances;
+
+  vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
+
+  bool distance1IsEnter = EnterDist == &Distance1;
+
+  vtkm::cont::ArrayHandle<FloatType> compactedDistance1;
+  vtkm::cont::Algorithm::CopyIf(Distance1, masks, compactedDistance1);
+  Distance1 = compactedDistance1;
+
+  vtkm::cont::ArrayHandle<FloatType> compactedDistance2;
+  vtkm::cont::Algorithm::CopyIf(Distance2, masks, compactedDistance2);
+  Distance2 = compactedDistance2;
+
+  vtkm::cont::ArrayHandle<vtkm::Int32> compactedExitFace;
+  vtkm::cont::Algorithm::CopyIf(ExitFace, masks, compactedExitFace);
+  ExitFace = compactedExitFace;
+
+  if (distance1IsEnter)
   {
-    tracer.RenderOnDevice(rays, device);
-    return true;
+    EnterDist = &Distance1;
+    ExitDist = &Distance2;
   }
-};
+  else
+  {
+    EnterDist = &Distance2;
+    ExitDist = &Distance1;
+  }
+}
+
+template <typename FloatType>
+void RayTracking<FloatType>::Init(const vtkm::Id size,
+                                  vtkm::cont::ArrayHandle<FloatType>& distances)
+{
+
+  ExitFace.Allocate(size);
+  Distance1.Allocate(size);
+  Distance2.Allocate(size);
+
+  CurrentDistance = distances;
+  //
+  // Set the initial Distances
+  //
+  vtkm::worklet::DispatcherMapField<CopyAndOffset<FloatType>> resetDistancesDispatcher(
+    CopyAndOffset<FloatType>(0.0f));
+  resetDistancesDispatcher.Invoke(distances, *EnterDist);
+
+  //
+  // Init the exit faces. This value is used to load the next cell
+  // base on the cell and face it left
+  //
+  vtkm::cont::ArrayHandleConstant<vtkm::Int32> negOne(-1, size);
+  vtkm::cont::Algorithm::Copy(negOne, ExitFace);
+
+  vtkm::cont::ArrayHandleConstant<FloatType> negOnef(-1.f, size);
+  vtkm::cont::Algorithm::Copy(negOnef, *ExitDist);
+}
+
+template <typename FloatType>
+void RayTracking<FloatType>::Swap()
+{
+  vtkm::cont::ArrayHandle<FloatType>* tmpPtr;
+  tmpPtr = EnterDist;
+  EnterDist = ExitDist;
+  ExitDist = tmpPtr;
+}
+
 } //namespace detail
 
-
-void ConnectivityTracer::Trace(Ray<vtkm::Float32>& rays)
-{
-  detail::RenderFunctor functor;
-  vtkm::cont::TryExecute(functor, *this, rays);
-}
-
-void ConnectivityTracer::Trace(Ray<vtkm::Float64>& rays)
-{
-  detail::RenderFunctor functor;
-  vtkm::cont::TryExecute(functor, *this, rays);
-}
-
 void ConnectivityTracer::Init()
 {
   //
@@ -220,9 +301,1179 @@ void ConnectivityTracer::LogTimers()
   logger->AddLogData("mesh_entry", LostRayTime);
 }
 
+template <typename FloatType>
+void ConnectivityTracer::PrintRayStatus(Ray<FloatType>& rays)
+{
+  vtkm::Id raysExited = RayOperations::GetStatusCount(rays, RAY_EXITED_MESH);
+  vtkm::Id raysActive = RayOperations::GetStatusCount(rays, RAY_ACTIVE);
+  vtkm::Id raysAbandoned = RayOperations::GetStatusCount(rays, RAY_ABANDONED);
+  vtkm::Id raysExitedDom = RayOperations::GetStatusCount(rays, RAY_EXITED_DOMAIN);
+  std::cout << "\r Ray Status " << std::setw(10) << std::left << " Lost " << std::setw(10)
+            << std::left << RaysLost << std::setw(10) << std::left << " Exited " << std::setw(10)
+            << std::left << raysExited << std::setw(10) << std::left << " Active " << std::setw(10)
+            << raysActive << std::setw(10) << std::left << " Abandoned " << std::setw(10)
+            << raysAbandoned << " Exited Domain " << std::setw(10) << std::left << raysExitedDom
+            << "\n";
+}
 
-template class detail::RayTracking<float>;
-template class detail::RayTracking<double>;
+//
+//  Advance Ray
+//      After a ray leaves the mesh, we need to check to see
+//      of the ray re-enters the mesh within this domain. This
+//      function moves the ray forward some offset to prevent
+//      "shadowing" and hitting the same exit point.
+//
+template <typename FloatType>
+class AdvanceRay : public vtkm::worklet::WorkletMapField
+{
+  FloatType Offset;
+
+public:
+  VTKM_CONT
+  AdvanceRay(const FloatType offset = 0.00001)
+    : Offset(offset)
+  {
+  }
+  using ControlSignature = void(FieldIn<>, FieldInOut<>);
+  using ExecutionSignature = void(_1, _2);
+
+  VTKM_EXEC inline void operator()(const vtkm::UInt8& status, FloatType& distance) const
+  {
+    if (status == RAY_EXITED_MESH)
+      distance += Offset;
+  }
+}; //class AdvanceRay
+
+class LocateCell : public vtkm::worklet::WorkletMapField
+{
+private:
+  CellIntersector<255> Intersector;
+
+public:
+  LocateCell() {}
+
+  using ControlSignature = void(FieldInOut<>,
+                                WholeArrayIn<>,
+                                FieldIn<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldIn<>,
+                                ExecObject meshConnectivity);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9);
+
+  template <typename FloatType, typename PointPortalType>
+  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
+                                   PointPortalType& vertices,
+                                   const vtkm::Vec<FloatType, 3>& dir,
+                                   FloatType& enterDistance,
+                                   FloatType& exitDistance,
+                                   vtkm::Int32& enterFace,
+                                   vtkm::UInt8& rayStatus,
+                                   const vtkm::Vec<FloatType, 3>& origin,
+                                   const MeshWrapper& meshConn) const
+  {
+    if (enterFace != -1 && rayStatus == RAY_ACTIVE)
+    {
+      currentCell = meshConn.GetConnectingCell(currentCell, enterFace);
+      if (currentCell == -1)
+        rayStatus = RAY_EXITED_MESH;
+      enterFace = -1;
+    }
+    //This ray is dead or exited the mesh and needs re-entry
+    if (rayStatus != RAY_ACTIVE)
+    {
+      return;
+    }
+    FloatType xpoints[8];
+    FloatType ypoints[8];
+    FloatType zpoints[8];
+    vtkm::Id cellConn[8];
+    FloatType distances[6];
+
+    const vtkm::Int32 numIndices = meshConn.GetCellIndices(cellConn, currentCell);
+    //load local cell data
+    for (int i = 0; i < numIndices; ++i)
+    {
+      BOUNDS_CHECK(vertices, cellConn[i]);
+      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
+      xpoints[i] = point[0];
+      ypoints[i] = point[1];
+      zpoints[i] = point[2];
+    }
+    const vtkm::UInt8 cellShape = meshConn.GetCellShape(currentCell);
+    Intersector.IntersectCell(xpoints, ypoints, zpoints, dir, origin, distances, cellShape);
+
+    CellTables tables;
+    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);
+    //vtkm::Int32 minFace = 6;
+    vtkm::Int32 maxFace = -1;
+
+    FloatType minDistance = static_cast<FloatType>(1e32);
+    FloatType maxDistance = static_cast<FloatType>(-1);
+    int hitCount = 0;
+    for (vtkm::Int32 i = 0; i < numFaces; ++i)
+    {
+      FloatType dist = distances[i];
+
+      if (dist != -1)
+      {
+        hitCount++;
+        if (dist < minDistance)
+        {
+          minDistance = dist;
+          //minFace = i;
+        }
+        if (dist > maxDistance)
+        {
+          maxDistance = dist;
+          maxFace = i;
+        }
+      }
+    }
+
+    if (maxDistance <= enterDistance || minDistance == maxDistance)
+    {
+      rayStatus = RAY_LOST;
+    }
+    else
+    {
+      enterDistance = minDistance;
+      exitDistance = maxDistance;
+      enterFace = maxFace;
+    }
+
+  } //operator
+};  //class LocateCell
+
+class RayBumper : public vtkm::worklet::WorkletMapField
+{
+private:
+  CellIntersector<255> Intersector;
+  const vtkm::UInt8 FailureStatus; // the status to assign ray if we fail to find the intersection
+public:
+  RayBumper(vtkm::UInt8 failureStatus = RAY_ABANDONED)
+    : FailureStatus(failureStatus)
+  {
+  }
+
+
+  using ControlSignature = void(FieldInOut<>,
+                                WholeArrayIn<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                FieldIn<>,
+                                FieldInOut<>,
+                                ExecObject meshConnectivity);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9);
+
+  template <typename FloatType, typename PointPortalType>
+  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
+                                   PointPortalType& vertices,
+                                   FloatType& enterDistance,
+                                   FloatType& exitDistance,
+                                   vtkm::Int32& enterFace,
+                                   vtkm::UInt8& rayStatus,
+                                   const vtkm::Vec<FloatType, 3>& origin,
+                                   vtkm::Vec<FloatType, 3>& rdir,
+                                   const MeshWrapper& meshConn) const
+  {
+    // We only process lost rays
+    if (rayStatus != RAY_LOST)
+    {
+      return;
+    }
+
+    FloatType xpoints[8];
+    FloatType ypoints[8];
+    FloatType zpoints[8];
+    vtkm::Id cellConn[8];
+    FloatType distances[6];
+
+    vtkm::Vec<FloatType, 3> centroid(0., 0., 0.);
+
+    const vtkm::Int32 numIndices = meshConn.GetCellIndices(cellConn, currentCell);
+    //load local cell data
+    for (int i = 0; i < numIndices; ++i)
+    {
+      BOUNDS_CHECK(vertices, cellConn[i]);
+      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
+      centroid = centroid + point;
+      xpoints[i] = point[0];
+      ypoints[i] = point[1];
+      zpoints[i] = point[2];
+    }
+
+    FloatType invNumIndices = static_cast<FloatType>(1.) / static_cast<FloatType>(numIndices);
+    centroid[0] = centroid[0] * invNumIndices;
+    centroid[1] = centroid[1] * invNumIndices;
+    centroid[2] = centroid[2] * invNumIndices;
+
+    vtkm::Vec<FloatType, 3> toCentroid = centroid - origin;
+    vtkm::Normalize(toCentroid);
+
+    vtkm::Vec<FloatType, 3> dir = rdir;
+    vtkm::Vec<FloatType, 3> bump = toCentroid - dir;
+    dir = dir + RAY_TUG_EPSILON * bump;
+
+    vtkm::Normalize(dir);
+    rdir = dir;
+
+    const vtkm::UInt8 cellShape = meshConn.GetCellShape(currentCell);
+    Intersector.IntersectCell(xpoints, ypoints, zpoints, rdir, origin, distances, cellShape);
+
+    CellTables tables;
+    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);
+
+    //vtkm::Int32 minFace = 6;
+    vtkm::Int32 maxFace = -1;
+    FloatType minDistance = static_cast<FloatType>(1e32);
+    FloatType maxDistance = static_cast<FloatType>(-1);
+    int hitCount = 0;
+    for (int i = 0; i < numFaces; ++i)
+    {
+      FloatType dist = distances[i];
+
+      if (dist != -1)
+      {
+        hitCount++;
+        if (dist < minDistance)
+        {
+          minDistance = dist;
+          //minFace = i;
+        }
+        if (dist >= maxDistance)
+        {
+          maxDistance = dist;
+          maxFace = i;
+        }
+      }
+    }
+    if (minDistance >= maxDistance)
+    {
+      rayStatus = FailureStatus;
+    }
+    else
+    {
+      enterDistance = minDistance;
+      exitDistance = maxDistance;
+      enterFace = maxFace;
+      rayStatus = RAY_ACTIVE; //re-activate ray
+    }
+
+  } //operator
+};  //class RayBumper
+
+class AddPathLengths : public vtkm::worklet::WorkletMapField
+{
+public:
+  VTKM_CONT
+  AddPathLengths() {}
+
+  using ControlSignature = void(FieldIn<RayStatusType>,            // ray status
+                                FieldIn<ScalarRenderingTypes>,     // cell enter distance
+                                FieldIn<ScalarRenderingTypes>,     // cell exit distance
+                                FieldInOut<ScalarRenderingTypes>); // ray absorption data
+
+  using ExecutionSignature = void(_1, _2, _3, _4);
+
+  template <typename FloatType>
+  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
+                                   const FloatType& enterDistance,
+                                   const FloatType& exitDistance,
+                                   FloatType& distance) const
+  {
+    if (rayStatus != RAY_ACTIVE)
+    {
+      return;
+    }
+
+    if (exitDistance <= enterDistance)
+    {
+      return;
+    }
+
+    FloatType segmentLength = exitDistance - enterDistance;
+    distance += segmentLength;
+  }
+};
+
+class Integrate : public vtkm::worklet::WorkletMapField
+{
+private:
+  const vtkm::Int32 NumBins;
+  const vtkm::Float32 UnitScalar;
+
+public:
+  VTKM_CONT
+  Integrate(const vtkm::Int32 numBins, const vtkm::Float32 unitScalar)
+    : NumBins(numBins)
+    , UnitScalar(unitScalar)
+  {
+  }
+
+  using ControlSignature = void(FieldIn<RayStatusType>,                // ray status
+                                FieldIn<ScalarRenderingTypes>,         // cell enter distance
+                                FieldIn<ScalarRenderingTypes>,         // cell exit distance
+                                FieldInOut<ScalarRenderingTypes>,      // current distance
+                                WholeArrayIn<ScalarRenderingTypes>,    // cell absorption data array
+                                WholeArrayInOut<ScalarRenderingTypes>, // ray absorption data
+                                FieldIn<IdType>);                      // current cell
+
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, WorkIndex);
+
+  template <typename FloatType, typename CellDataPortalType, typename RayDataPortalType>
+  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
+                                   const FloatType& enterDistance,
+                                   const FloatType& exitDistance,
+                                   FloatType& currentDistance,
+                                   const CellDataPortalType& cellData,
+                                   RayDataPortalType& energyBins,
+                                   const vtkm::Id& currentCell,
+                                   const vtkm::Id& rayIndex) const
+  {
+    if (rayStatus != RAY_ACTIVE)
+    {
+      return;
+    }
+    if (exitDistance <= enterDistance)
+    {
+      return;
+    }
+
+    FloatType segmentLength = exitDistance - enterDistance;
+
+    vtkm::Id rayOffset = NumBins * rayIndex;
+    vtkm::Id cellOffset = NumBins * currentCell;
+    for (vtkm::Int32 i = 0; i < NumBins; ++i)
+    {
+      BOUNDS_CHECK(cellData, cellOffset + i);
+      FloatType absorb = static_cast<FloatType>(cellData.Get(cellOffset + i));
+      absorb *= UnitScalar;
+      absorb = vtkm::Exp(-absorb * segmentLength);
+      BOUNDS_CHECK(energyBins, rayOffset + i);
+      FloatType intensity = static_cast<FloatType>(energyBins.Get(rayOffset + i));
+      energyBins.Set(rayOffset + i, intensity * absorb);
+    }
+    currentDistance = exitDistance;
+  }
+};
+
+class IntegrateEmission : public vtkm::worklet::WorkletMapField
+{
+private:
+  const vtkm::Int32 NumBins;
+  const vtkm::Float32 UnitScalar;
+  bool DivideEmisByAbsorb;
+
+public:
+  VTKM_CONT
+  IntegrateEmission(const vtkm::Int32 numBins,
+                    const vtkm::Float32 unitScalar,
+                    const bool divideEmisByAbsorb)
+    : NumBins(numBins)
+    , UnitScalar(unitScalar)
+    , DivideEmisByAbsorb(divideEmisByAbsorb)
+  {
+  }
+
+  using ControlSignature = void(FieldIn<>,                          // ray status
+                                FieldIn<>,                          // cell enter distance
+                                FieldIn<>,                          // cell exit distance
+                                FieldInOut<>,                       // current distance
+                                WholeArrayIn<ScalarRenderingTypes>, // cell absorption data array
+                                WholeArrayIn<ScalarRenderingTypes>, // cell emission data array
+                                WholeArrayInOut<>,                  // ray absorption data
+                                WholeArrayInOut<>,                  // ray emission data
+                                FieldIn<>);                         // current cell
+
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9, WorkIndex);
+
+  template <typename FloatType,
+            typename CellAbsPortalType,
+            typename CellEmisPortalType,
+            typename RayDataPortalType>
+  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
+                                   const FloatType& enterDistance,
+                                   const FloatType& exitDistance,
+                                   FloatType& currentDistance,
+                                   const CellAbsPortalType& absorptionData,
+                                   const CellEmisPortalType& emissionData,
+                                   RayDataPortalType& absorptionBins,
+                                   RayDataPortalType& emissionBins,
+                                   const vtkm::Id& currentCell,
+                                   const vtkm::Id& rayIndex) const
+  {
+    if (rayStatus != RAY_ACTIVE)
+    {
+      return;
+    }
+    if (exitDistance <= enterDistance)
+    {
+      return;
+    }
+
+    FloatType segmentLength = exitDistance - enterDistance;
+
+    vtkm::Id rayOffset = NumBins * rayIndex;
+    vtkm::Id cellOffset = NumBins * currentCell;
+    for (vtkm::Int32 i = 0; i < NumBins; ++i)
+    {
+      BOUNDS_CHECK(absorptionData, cellOffset + i);
+      FloatType absorb = static_cast<FloatType>(absorptionData.Get(cellOffset + i));
+      BOUNDS_CHECK(emissionData, cellOffset + i);
+      FloatType emission = static_cast<FloatType>(emissionData.Get(cellOffset + i));
+
+      absorb *= UnitScalar;
+      emission *= UnitScalar;
+
+      if (DivideEmisByAbsorb)
+      {
+        emission /= absorb;
+      }
+
+      FloatType tmp = vtkm::Exp(-absorb * segmentLength);
+      BOUNDS_CHECK(absorptionBins, rayOffset + i);
+
+      //
+      // Traditionally, we would only keep track of a single intensity value per ray
+      // per bin and we would integrate from the beginning to end of the ray. In a
+      // distributed memory setting, we would move cell data around so that the
+      // entire ray could be traced, but in situ, moving that much cell data around
+      // could blow memory. Here we are keeping track of two values. Total absorption
+      // through this contiguous segment of the mesh, and the amount of emitted energy
+      // that makes it out of this mesh segment. If this is really run on a single node,
+      // we can get the final energy value by multiplying the background intensity by
+      // the total absorption of the mesh segment and add in the amount of emitted
+      // energy that escapes.
+      //
+      FloatType absorbIntensity = static_cast<FloatType>(absorptionBins.Get(rayOffset + i));
+      FloatType emissionIntensity = static_cast<FloatType>(emissionBins.Get(rayOffset + i));
+
+      absorptionBins.Set(rayOffset + i, absorbIntensity * tmp);
+
+      emissionIntensity = emissionIntensity * tmp + emission * (1.f - tmp);
+
+      BOUNDS_CHECK(emissionBins, rayOffset + i);
+      emissionBins.Set(rayOffset + i, emissionIntensity);
+    }
+    currentDistance = exitDistance;
+  }
+};
+//
+//  IdentifyMissedRay is a debugging routine that detects
+//  rays that fail to have any value because of a external
+//  intersection and cell intersection mismatch
+//
+//
+class IdentifyMissedRay : public vtkm::worklet::WorkletMapField
+{
+public:
+  vtkm::Id Width;
+  vtkm::Id Height;
+  vtkm::Vec<vtkm::Float32, 4> BGColor;
+  IdentifyMissedRay(const vtkm::Id width,
+                    const vtkm::Id height,
+                    vtkm::Vec<vtkm::Float32, 4> bgcolor)
+    : Width(width)
+    , Height(height)
+    , BGColor(bgcolor)
+  {
+  }
+  using ControlSignature = void(FieldIn<>, WholeArrayIn<>);
+  using ExecutionSignature = void(_1, _2);
+
+
+  VTKM_EXEC inline bool IsBGColor(const vtkm::Vec<vtkm::Float32, 4> color) const
+  {
+    bool isBG = false;
+
+    if (color[0] == BGColor[0] && color[1] == BGColor[1] && color[2] == BGColor[2] &&
+        color[3] == BGColor[3])
+      isBG = true;
+    return isBG;
+  }
+
+  template <typename ColorBufferType>
+  VTKM_EXEC inline void operator()(const vtkm::Id& pixelId, ColorBufferType& buffer) const
+  {
+    vtkm::Id x = pixelId % Width;
+    vtkm::Id y = pixelId / Width;
+
+    // Conservative check, we only want to check pixels in the middle
+    if (x <= 0 || y <= 0)
+      return;
+    if (x >= Width - 1 || y >= Height - 1)
+      return;
+    vtkm::Vec<vtkm::Float32, 4> pixel;
+    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 0));
+    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 1));
+    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 2));
+    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 3));
+    if (!IsBGColor(pixel))
+      return;
+    vtkm::Id p0 = (y)*Width + (x + 1);
+    vtkm::Id p1 = (y)*Width + (x - 1);
+    vtkm::Id p2 = (y + 1) * Width + (x);
+    vtkm::Id p3 = (y - 1) * Width + (x);
+    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 0));
+    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 1));
+    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 2));
+    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 3));
+    if (IsBGColor(pixel))
+      return;
+    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 0));
+    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 1));
+    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 2));
+    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 3));
+    if (IsBGColor(pixel))
+      return;
+    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 0));
+    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 1));
+    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 2));
+    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 3));
+    if (IsBGColor(pixel))
+      return;
+    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 0));
+    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 1));
+    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 2));
+    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 3));
+    if (IsBGColor(pixel))
+      return;
+
+    printf("Possible error ray missed ray %d\n", (int)pixelId);
+  }
+};
+
+template <typename FloatType>
+class SampleCellAssocCells : public vtkm::worklet::WorkletMapField
+{
+private:
+  CellSampler<255> Sampler;
+  FloatType SampleDistance;
+  FloatType MinScalar;
+  FloatType InvDeltaScalar;
+
+public:
+  SampleCellAssocCells(const FloatType& sampleDistance,
+                       const FloatType& minScalar,
+                       const FloatType& maxScalar)
+    : SampleDistance(sampleDistance)
+    , MinScalar(minScalar)
+  {
+    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
+  }
+
+
+  using ControlSignature = void(FieldIn<>,
+                                WholeArrayIn<ScalarRenderingTypes>,
+                                FieldIn<>,
+                                FieldIn<>,
+                                FieldInOut<>,
+                                FieldInOut<>,
+                                WholeArrayIn<>,
+                                WholeArrayInOut<>);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, WorkIndex);
+
+  template <typename ScalarPortalType, typename ColorMapType, typename FrameBufferType>
+  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
+                                   ScalarPortalType& scalarPortal,
+                                   const FloatType& enterDistance,
+                                   const FloatType& exitDistance,
+                                   FloatType& currentDistance,
+                                   vtkm::UInt8& rayStatus,
+                                   const ColorMapType& colorMap,
+                                   FrameBufferType& frameBuffer,
+                                   const vtkm::Id& pixelIndex) const
+  {
+
+    if (rayStatus != RAY_ACTIVE)
+      return;
+
+    vtkm::Vec<vtkm::Float32, 4> color;
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 0);
+    color[0] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 0));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 1);
+    color[1] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 1));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 2);
+    color[2] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 2));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 3);
+    color[3] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 3));
+
+    vtkm::Float32 scalar;
+    BOUNDS_CHECK(scalarPortal, currentCell);
+    scalar = vtkm::Float32(scalarPortal.Get(currentCell));
+    //
+    // There can be mismatches in the initial enter distance and the current distance
+    // due to lost rays at cell borders. For now,
+    // we will just advance the current position to the enter distance, since otherwise,
+    // the pixel would never be sampled.
+    //
+    if (currentDistance < enterDistance)
+      currentDistance = enterDistance;
+
+    const vtkm::Id colorMapSize = colorMap.GetNumberOfValues();
+    vtkm::Float32 lerpedScalar;
+    lerpedScalar = static_cast<vtkm::Float32>((scalar - MinScalar) * InvDeltaScalar);
+    vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(colorMapSize));
+    if (colorIndex < 0)
+      colorIndex = 0;
+    if (colorIndex >= colorMapSize)
+      colorIndex = colorMapSize - 1;
+    BOUNDS_CHECK(colorMap, colorIndex);
+    vtkm::Vec<vtkm::Float32, 4> sampleColor = colorMap.Get(colorIndex);
+
+    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
+    {
+      //composite
+      vtkm::Float32 alpha = sampleColor[3] * (1.f - color[3]);
+      color[0] = color[0] + sampleColor[0] * alpha;
+      color[1] = color[1] + sampleColor[1] * alpha;
+      color[2] = color[2] + sampleColor[2] * alpha;
+      color[3] = alpha + color[3];
+
+      if (color[3] > 1.)
+      {
+        rayStatus = RAY_TERMINATED;
+        break;
+      }
+      currentDistance += SampleDistance;
+    }
+
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 0);
+    frameBuffer.Set(pixelIndex * 4 + 0, color[0]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 1);
+    frameBuffer.Set(pixelIndex * 4 + 1, color[1]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 2);
+    frameBuffer.Set(pixelIndex * 4 + 2, color[2]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 3);
+    frameBuffer.Set(pixelIndex * 4 + 3, color[3]);
+  }
+}; //class Sample cell
+
+template <typename FloatType>
+class SampleCellAssocPoints : public vtkm::worklet::WorkletMapField
+{
+private:
+  CellSampler<255> Sampler;
+  FloatType SampleDistance;
+  FloatType MinScalar;
+  FloatType InvDeltaScalar;
+
+public:
+  SampleCellAssocPoints(const FloatType& sampleDistance,
+                        const FloatType& minScalar,
+                        const FloatType& maxScalar)
+    : SampleDistance(sampleDistance)
+    , MinScalar(minScalar)
+  {
+    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
+  }
+
+
+  using ControlSignature = void(FieldIn<>,
+                                WholeArrayIn<Vec3>,
+                                WholeArrayIn<ScalarRenderingTypes>,
+                                FieldIn<>,
+                                FieldIn<>,
+                                FieldInOut<>,
+                                FieldIn<>,
+                                FieldInOut<>,
+                                FieldIn<>,
+                                ExecObject meshConnectivity,
+                                WholeArrayIn<>,
+                                WholeArrayInOut<>);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, WorkIndex, _9, _10, _11, _12);
+
+  template <typename PointPortalType,
+            typename ScalarPortalType,
+            typename ColorMapType,
+            typename FrameBufferType>
+  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
+                                   PointPortalType& vertices,
+                                   ScalarPortalType& scalarPortal,
+                                   const FloatType& enterDistance,
+                                   const FloatType& exitDistance,
+                                   FloatType& currentDistance,
+                                   const vtkm::Vec<vtkm::Float32, 3>& dir,
+                                   vtkm::UInt8& rayStatus,
+                                   const vtkm::Id& pixelIndex,
+                                   const vtkm::Vec<FloatType, 3>& origin,
+                                   MeshWrapper& meshConn,
+                                   const ColorMapType& colorMap,
+                                   FrameBufferType& frameBuffer) const
+  {
+
+    if (rayStatus != RAY_ACTIVE)
+      return;
+
+    vtkm::Vec<vtkm::Float32, 4> color;
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 0);
+    color[0] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 0));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 1);
+    color[1] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 1));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 2);
+    color[2] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 2));
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 3);
+    color[3] = static_cast<vtkm::Float32>(frameBuffer.Get(pixelIndex * 4 + 3));
+
+    if (color[3] >= 1.f)
+    {
+      rayStatus = RAY_TERMINATED;
+      return;
+    }
+    vtkm::Vec<vtkm::Float32, 8> scalars;
+    vtkm::Vec<vtkm::Vec<FloatType, 3>, 8> points;
+    // silence "may" be uninitialized warning
+    for (vtkm::Int32 i = 0; i < 8; ++i)
+    {
+      scalars[i] = 0.f;
+      points[i] = vtkm::Vec<FloatType, 3>(0.f, 0.f, 0.f);
+    }
+    //load local scalar cell data
+    vtkm::Id cellConn[8];
+    const vtkm::Int32 numIndices = meshConn.GetCellIndices(cellConn, currentCell);
+    for (int i = 0; i < numIndices; ++i)
+    {
+      BOUNDS_CHECK(scalarPortal, cellConn[i]);
+      scalars[i] = static_cast<vtkm::Float32>(scalarPortal.Get(cellConn[i]));
+      BOUNDS_CHECK(vertices, cellConn[i]);
+      points[i] = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
+    }
+    //
+    // There can be mismatches in the initial enter distance and the current distance
+    // due to lost rays at cell borders. For now,
+    // we will just advance the current position to the enter distance, since otherwise,
+    // the pixel would never be sampled.
+    //
+    if (currentDistance < enterDistance)
+    {
+      currentDistance = enterDistance;
+    }
+
+    const vtkm::Id colorMapSize = colorMap.GetNumberOfValues();
+    const vtkm::Int32 cellShape = meshConn.GetCellShape(currentCell);
+
+    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
+    {
+      vtkm::Vec<FloatType, 3> sampleLoc = origin + currentDistance * dir;
+      vtkm::Float32 lerpedScalar;
+      bool validSample =
+        Sampler.SampleCell(points, scalars, sampleLoc, lerpedScalar, *this, cellShape);
+      if (!validSample)
+      {
+        //
+        // There is a slight mismatch between intersections and parametric coordinates
+        // which results in a invalid sample very close to the cell edge. Just throw
+        // this sample away, and move to the next sample.
+        //
+
+        //There should be a sample here, so offset and try again.
+
+        currentDistance += 0.00001f;
+        continue;
+      }
+      lerpedScalar = static_cast<vtkm::Float32>((lerpedScalar - MinScalar) * InvDeltaScalar);
+      vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(colorMapSize));
+
+      colorIndex = vtkm::Min(vtkm::Max(colorIndex, vtkm::Id(0)), colorMapSize - 1);
+      BOUNDS_CHECK(colorMap, colorIndex);
+      vtkm::Vec<vtkm::Float32, 4> sampleColor = colorMap.Get(colorIndex);
+      //composite
+      sampleColor[3] *= (1.f - color[3]);
+      color[0] = color[0] + sampleColor[0] * sampleColor[3];
+      color[1] = color[1] + sampleColor[1] * sampleColor[3];
+      color[2] = color[2] + sampleColor[2] * sampleColor[3];
+      color[3] = sampleColor[3] + color[3];
+
+      if (color[3] >= 1.0)
+      {
+        rayStatus = RAY_TERMINATED;
+        break;
+      }
+      currentDistance += SampleDistance;
+    }
+
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 0);
+    frameBuffer.Set(pixelIndex * 4 + 0, color[0]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 1);
+    frameBuffer.Set(pixelIndex * 4 + 1, color[1]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 2);
+    frameBuffer.Set(pixelIndex * 4 + 2, color[2]);
+    BOUNDS_CHECK(frameBuffer, pixelIndex * 4 + 3);
+    frameBuffer.Set(pixelIndex * 4 + 3, color[3]);
+  }
+}; //class Sample cell
+
+template <typename FloatType>
+void ConnectivityTracer::IntersectCell(Ray<FloatType>& rays,
+                                       detail::RayTracking<FloatType>& tracker)
+{
+  vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
+  vtkm::worklet::DispatcherMapField<LocateCell> locateDispatch;
+  locateDispatch.Invoke(rays.HitIdx,
+                        this->Coords,
+                        rays.Dir,
+                        *(tracker.EnterDist),
+                        *(tracker.ExitDist),
+                        tracker.ExitFace,
+                        rays.Status,
+                        rays.Origin,
+                        MeshContainer);
+
+  if (this->CountRayStatus)
+    RaysLost = RayOperations::GetStatusCount(rays, RAY_LOST);
+  this->IntersectTime += timer.GetElapsedTime();
+}
+
+template <typename FloatType>
+void ConnectivityTracer::AccumulatePathLengths(Ray<FloatType>& rays,
+                                               detail::RayTracking<FloatType>& tracker)
+{
+  vtkm::worklet::DispatcherMapField<AddPathLengths> dispatcher;
+  dispatcher.Invoke(
+    rays.Status, *(tracker.EnterDist), *(tracker.ExitDist), rays.GetBuffer("path_lengths").Buffer);
+}
+
+template <typename FloatType>
+void ConnectivityTracer::FindLostRays(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker)
+{
+  vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
+
+  vtkm::worklet::DispatcherMapField<RayBumper> bumpDispatch;
+  bumpDispatch.Invoke(rays.HitIdx,
+                      this->Coords,
+                      *(tracker.EnterDist),
+                      *(tracker.ExitDist),
+                      tracker.ExitFace,
+                      rays.Status,
+                      rays.Origin,
+                      rays.Dir,
+                      MeshContainer);
+
+  this->LostRayTime += timer.GetElapsedTime();
+}
+
+template <typename FloatType>
+void ConnectivityTracer::SampleCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker)
+{
+  using SampleP = SampleCellAssocPoints<FloatType>;
+  using SampleC = SampleCellAssocCells<FloatType>;
+  vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
+
+  VTKM_ASSERT(rays.Buffers.at(0).GetNumChannels() == 4);
+
+  if (FieldAssocPoints)
+  {
+    vtkm::worklet::DispatcherMapField<SampleP> dispatcher(
+      SampleP(this->SampleDistance,
+              vtkm::Float32(this->ScalarBounds.Min),
+              vtkm::Float32(this->ScalarBounds.Max)));
+    dispatcher.Invoke(rays.HitIdx,
+                      this->Coords,
+                      this->ScalarField.GetData(),
+                      *(tracker.EnterDist),
+                      *(tracker.ExitDist),
+                      tracker.CurrentDistance,
+                      rays.Dir,
+                      rays.Status,
+                      rays.Origin,
+                      MeshContainer,
+                      this->ColorMap,
+                      rays.Buffers.at(0).Buffer);
+  }
+  else
+  {
+    vtkm::worklet::DispatcherMapField<SampleC> dispatcher(
+      SampleC(this->SampleDistance,
+              vtkm::Float32(this->ScalarBounds.Min),
+              vtkm::Float32(this->ScalarBounds.Max)));
+
+    dispatcher.Invoke(rays.HitIdx,
+                      this->ScalarField.GetData(),
+                      *(tracker.EnterDist),
+                      *(tracker.ExitDist),
+                      tracker.CurrentDistance,
+                      rays.Status,
+                      this->ColorMap,
+                      rays.Buffers.at(0).Buffer);
+  }
+
+  this->SampleTime += timer.GetElapsedTime();
+}
+
+template <typename FloatType>
+void ConnectivityTracer::IntegrateCells(Ray<FloatType>& rays,
+                                        detail::RayTracking<FloatType>& tracker)
+{
+  vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
+  if (HasEmission)
+  {
+    bool divideEmisByAbsorp = false;
+    vtkm::cont::ArrayHandle<FloatType> absorp = rays.Buffers.at(0).Buffer;
+    vtkm::cont::ArrayHandle<FloatType> emission = rays.GetBuffer("emission").Buffer;
+    vtkm::worklet::DispatcherMapField<IntegrateEmission> dispatcher(
+      IntegrateEmission(rays.Buffers.at(0).GetNumChannels(), UnitScalar, divideEmisByAbsorp));
+    dispatcher.Invoke(rays.Status,
+                      *(tracker.EnterDist),
+                      *(tracker.ExitDist),
+                      rays.Distance,
+                      this->ScalarField.GetData(),
+                      this->EmissionField.GetData(),
+                      absorp,
+                      emission,
+                      rays.HitIdx);
+  }
+  else
+  {
+    vtkm::worklet::DispatcherMapField<Integrate> dispatcher(
+      Integrate(rays.Buffers.at(0).GetNumChannels(), UnitScalar));
+    dispatcher.Invoke(rays.Status,
+                      *(tracker.EnterDist),
+                      *(tracker.ExitDist),
+                      rays.Distance,
+                      this->ScalarField.GetData(),
+                      rays.Buffers.at(0).Buffer,
+                      rays.HitIdx);
+  }
+
+  IntegrateTime += timer.GetElapsedTime();
+}
+
+// template <typename FloatType>
+// void ConnectivityTracer<CellType>::PrintDebugRay(Ray<FloatType>& rays, vtkm::Id rayId)
+// {
+//   vtkm::Id index = -1;
+//   for (vtkm::Id i = 0; i < rays.NumRays; ++i)
+//   {
+//     if (rays.PixelIdx.GetPortalControl().Get(i) == rayId)
+//     {
+//       index = i;
+//       break;
+//     }
+//   }
+//   if (index == -1)
+//   {
+//     return;
+//   }
+
+//   std::cout << "++++++++RAY " << rayId << "++++++++\n";
+//   std::cout << "Status: " << (int)rays.Status.GetPortalControl().Get(index) << "\n";
+//   std::cout << "HitIndex: " << rays.HitIdx.GetPortalControl().Get(index) << "\n";
+//   std::cout << "Dist " << rays.Distance.GetPortalControl().Get(index) << "\n";
+//   std::cout << "MinDist " << rays.MinDistance.GetPortalControl().Get(index) << "\n";
+//   std::cout << "Origin " << rays.Origin.GetPortalConstControl().Get(index) << "\n";
+//   std::cout << "Dir " << rays.Dir.GetPortalConstControl().Get(index) << "\n";
+//   std::cout << "+++++++++++++++++++++++++\n";
+// }
+
+template <typename FloatType>
+void ConnectivityTracer::OffsetMinDistances(Ray<FloatType>& rays)
+{
+  vtkm::worklet::DispatcherMapField<AdvanceRay<FloatType>> dispatcher(
+    AdvanceRay<FloatType>(FloatType(0.001)));
+  dispatcher.Invoke(rays.Status, rays.MinDistance);
+}
+
+template <typename FloatType>
+void ConnectivityTracer::FindMeshEntry(Ray<FloatType>& rays)
+{
+  vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> entryTimer;
+  //
+  // if ray misses the external face it will be marked RAY_EXITED_MESH
+  //
+  MeshContainer->FindEntry(rays);
+  MeshEntryTime += entryTimer.GetElapsedTime();
+}
+
+template <typename FloatType>
+void ConnectivityTracer::IntegrateMeshSegment(Ray<FloatType>& rays)
+{
+  this->Init(); // sets sample distance
+  detail::RayTracking<FloatType> rayTracker;
+  rayTracker.Init(rays.NumRays, rays.Distance);
+
+  bool hasPathLengths = rays.HasBuffer("path_lengths");
+
+  while (RayOperations::RaysInMesh(rays))
+  {
+    //
+    // Rays the leave the mesh will be marked as RAYEXITED_MESH
+    this->IntersectCell(rays, rayTracker);
+    //
+    // If the ray was lost due to precision issues, we find it.
+    // If it is marked RAY_ABANDONED, then something went wrong.
+    //
+    this->FindLostRays(rays, rayTracker);
+    //
+    // integrate along the ray
+    //
+    if (this->Integrator == Volume)
+      this->SampleCells(rays, rayTracker);
+    else
+      this->IntegrateCells(rays, rayTracker);
+
+    if (hasPathLengths)
+    {
+      this->AccumulatePathLengths(rays, rayTracker);
+    }
+    //swap enter and exit distances
+    rayTracker.Swap();
+    if (this->CountRayStatus)
+      this->PrintRayStatus(rays);
+  } //for
+}
+
+template <typename FloatType>
+void ConnectivityTracer::FullTrace(Ray<FloatType>& rays)
+{
+
+  this->RaysLost = 0;
+  RayOperations::ResetStatus(rays, RAY_EXITED_MESH);
+
+  if (this->CountRayStatus)
+  {
+    this->PrintRayStatus(rays);
+  }
+
+  bool cullMissedRays = true;
+  bool workRemaining = true;
+
+  do
+  {
+    FindMeshEntry(rays);
+
+    if (cullMissedRays)
+    {
+      vtkm::cont::ArrayHandle<UInt8> activeRays;
+      activeRays = RayOperations::CompactActiveRays(rays);
+      cullMissedRays = false;
+    }
+
+    IntegrateMeshSegment(rays);
+
+    workRemaining = RayOperations::RaysProcessed(rays) != rays.NumRays;
+    //
+    // Ensure that we move the current distance forward some
+    // epsilon so we don't re-enter the cell we just left.
+    //
+    if (workRemaining)
+    {
+      RayOperations::CopyDistancesToMin(rays);
+      this->OffsetMinDistances(rays);
+    }
+  } while (workRemaining);
+}
+
+template <typename FloatType>
+std::vector<PartialComposite<FloatType>> ConnectivityTracer::PartialTrace(Ray<FloatType>& rays)
+{
+
+  bool hasPathLengths = rays.HasBuffer("path_lengths");
+  this->RaysLost = 0;
+  RayOperations::ResetStatus(rays, RAY_EXITED_MESH);
+
+  std::vector<PartialComposite<FloatType>> partials;
+
+  if (this->CountRayStatus)
+  {
+    this->PrintRayStatus(rays);
+  }
+
+  bool workRemaining = true;
+
+  do
+  {
+    FindMeshEntry(rays);
+
+    vtkm::cont::ArrayHandle<UInt8> activeRays;
+    activeRays = RayOperations::CompactActiveRays(rays);
+
+    if (rays.NumRays == 0)
+      break;
+
+    IntegrateMeshSegment(rays);
+
+    PartialComposite<FloatType> partial;
+    partial.Buffer = rays.Buffers.at(0).Copy();
+    vtkm::cont::Algorithm::Copy(rays.Distance, partial.Distances);
+    vtkm::cont::Algorithm::Copy(rays.PixelIdx, partial.PixelIds);
+
+    if (HasEmission && this->Integrator == Energy)
+    {
+      partial.Intensities = rays.GetBuffer("emission").Copy();
+    }
+    if (hasPathLengths)
+    {
+      partial.PathLengths = rays.GetBuffer("path_lengths").Copy().Buffer;
+    }
+    partials.push_back(partial);
+
+    // reset buffers
+    if (this->Integrator == Volume)
+    {
+      vtkm::cont::ArrayHandle<FloatType> signature;
+      signature.Allocate(4);
+      signature.GetPortalControl().Set(0, 0.f);
+      signature.GetPortalControl().Set(1, 0.f);
+      signature.GetPortalControl().Set(2, 0.f);
+      signature.GetPortalControl().Set(3, 0.f);
+      rays.Buffers.at(0).InitChannels(signature);
+    }
+    else
+    {
+      rays.Buffers.at(0).InitConst(1.f);
+      if (HasEmission)
+      {
+        rays.GetBuffer("emission").InitConst(0.f);
+      }
+      if (hasPathLengths)
+      {
+        rays.GetBuffer("path_lengths").InitConst(0.f);
+      }
+    }
+
+    workRemaining = RayOperations::RaysProcessed(rays) != rays.NumRays;
+    //
+    // Ensure that we move the current distance forward some
+    // epsilon so we don't re-enter the cell we just left.
+    //
+    if (workRemaining)
+    {
+      RayOperations::CopyDistancesToMin(rays);
+      this->OffsetMinDistances(rays);
+    }
+  } while (workRemaining);
+
+  return partials;
+}
+
+template class detail::RayTracking<vtkm::Float32>;
+template class detail::RayTracking<vtkm::Float64>;
+
+template struct PartialComposite<vtkm::Float32>;
+template struct PartialComposite<vtkm::Float64>;
+
+template void ConnectivityTracer::FullTrace<vtkm::Float32>(Ray<vtkm::Float32>& rays);
+
+template std::vector<PartialComposite<vtkm::Float32>>
+ConnectivityTracer::PartialTrace<vtkm::Float32>(Ray<vtkm::Float32>& rays);
+
+template void ConnectivityTracer::IntegrateMeshSegment<vtkm::Float32>(Ray<vtkm::Float32>& rays);
+
+template void ConnectivityTracer::FindMeshEntry<vtkm::Float32>(Ray<vtkm::Float32>& rays);
+
+template void ConnectivityTracer::FullTrace<vtkm::Float64>(Ray<vtkm::Float64>& rays);
+
+template std::vector<PartialComposite<vtkm::Float64>>
+ConnectivityTracer::PartialTrace<vtkm::Float64>(Ray<vtkm::Float64>& rays);
+
+template void ConnectivityTracer::IntegrateMeshSegment<vtkm::Float64>(Ray<vtkm::Float64>& rays);
+
+template void ConnectivityTracer::FindMeshEntry<vtkm::Float64>(Ray<vtkm::Float64>& rays);
 }
 }
 } // namespace vtkm::rendering::raytracing
diff --git a/vtkm/rendering/raytracing/ConnectivityTracer.h b/vtkm/rendering/raytracing/ConnectivityTracer.h
index 831e0395e..9a0c6c4bb 100644
--- a/vtkm/rendering/raytracing/ConnectivityTracer.h
+++ b/vtkm/rendering/raytracing/ConnectivityTracer.h
@@ -22,18 +22,11 @@
 
 #include <vtkm/rendering/vtkm_rendering_export.h>
 
-#include <iomanip>
 #include <vtkm/cont/ArrayHandle.h>
 
 #include <vtkm/rendering/raytracing/MeshConnectivityContainers.h>
+#include <vtkm/rendering/raytracing/PartialComposite.h>
 
-#ifndef CELL_SHAPE_ZOO
-#define CELL_SHAPE_ZOO 255
-#endif
-
-#ifndef CELL_SHAPE_STRUCTURED
-#define CELL_SHAPE_STRUCTURED 254
-#endif
 
 namespace vtkm
 {
@@ -68,26 +61,29 @@ public:
     ExitDist = &Distance2;
   }
 
-  template <typename Device>
   void Compact(vtkm::cont::ArrayHandle<FloatType>& compactedDistances,
-               vtkm::cont::ArrayHandle<UInt8>& masks,
-               Device);
+               vtkm::cont::ArrayHandle<UInt8>& masks);
 
-  template <typename Device>
-  void Init(const vtkm::Id size, vtkm::cont::ArrayHandle<FloatType>& distances, Device);
+  void Init(const vtkm::Id size, vtkm::cont::ArrayHandle<FloatType>& distances);
 
   void Swap();
 };
 
 } //namespace detail
 
-
-class ConnectivityTracer
+/**
+ * \brief ConnectivityTracer is volumetric ray tracer for unstructured
+ *        grids. Capabilities include volume rendering and integrating
+ *        absorption and emission of N energy groups for simulated
+ *        radiograhy.
+ */
+class VTKM_RENDERING_EXPORT ConnectivityTracer
 {
 public:
   ConnectivityTracer()
     : MeshContainer(nullptr)
     , CountRayStatus(false)
+    , UnitScalar(1.f)
   {
   }
 
@@ -120,57 +116,75 @@ public:
   void SetSampleDistance(const vtkm::Float32& distance);
   void SetColorMap(const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>& colorMap);
 
-  void Trace(Ray<vtkm::Float32>& rays);
-
-  void Trace(Ray<vtkm::Float64>& rays);
-
   MeshConnContainer* GetMeshContainer() { return MeshContainer; }
 
   void Init();
 
   void SetDebugOn(bool on) { CountRayStatus = on; }
 
+  void SetUnitScalar(const vtkm::Float32 unitScalar) { UnitScalar = unitScalar; }
+
+
   vtkm::Id GetNumberOfMeshCells() const;
 
   void ResetTimers();
   void LogTimers();
 
+  ///
+  /// Traces rays fully through the mesh. Rays can exit and re-enter
+  /// multiple times before leaving the domain. This is fast path for
+  /// structured meshs or meshes that are not interlocking.
+  /// Note: rays will be compacted
+  ///
+  template <typename FloatType>
+  void FullTrace(Ray<FloatType>& rays);
+
+  ///
+  /// Integrates rays through the mesh. If rays leave the mesh and
+  /// re-enter, then those become two separate partial composites.
+  /// This is need to support domain decompositions that are like
+  /// puzzle pieces. Note: rays will be compacted
+  ///
+  template <typename FloatType>
+  std::vector<PartialComposite<FloatType>> PartialTrace(Ray<FloatType>& rays);
+
+  ///
+  /// Integrates the active rays though the mesh until all rays
+  /// have exited.
+  ///  Precondition: rays.HitIdx is set to a valid mesh cell
+  ///
+  template <typename FloatType>
+  void IntegrateMeshSegment(Ray<FloatType>& rays);
+
+  ///
+  /// Find the entry point in the mesh
+  ///
+  template <typename FloatType>
+  void FindMeshEntry(Ray<FloatType>& rays);
+
 private:
-  friend struct detail::RenderFunctor;
-  template <typename FloatType, typename Device>
-  void IntersectCell(Ray<FloatType>& rays,
-                     detail::RayTracking<FloatType>& tracker,
-                     const MeshConnectivityBase* meshConn,
-                     Device);
+  template <typename FloatType>
+  void IntersectCell(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
 
-  template <typename FloatType, typename Device>
-  void AccumulatePathLengths(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker, Device);
+  template <typename FloatType>
+  void AccumulatePathLengths(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
 
-  template <typename FloatType, typename Device>
-  void FindLostRays(Ray<FloatType>& rays,
-                    detail::RayTracking<FloatType>& tracker,
-                    const MeshConnectivityBase* meshConn,
-                    Device);
+  template <typename FloatType>
+  void FindLostRays(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
 
-  template <typename FloatType, typename Device>
-  void SampleCells(Ray<FloatType>& rays,
-                   detail::RayTracking<FloatType>& tracker,
-                   const MeshConnectivityBase* meshConn,
-                   Device);
+  template <typename FloatType>
+  void SampleCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
 
-  template <typename FloatType, typename Device>
-  void IntegrateCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker, Device);
+  template <typename FloatType>
+  void IntegrateCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
 
-  template <typename FloatType, typename Device>
-  void OffsetMinDistances(Ray<FloatType>& rays, Device);
+  template <typename FloatType>
+  void OffsetMinDistances(Ray<FloatType>& rays);
 
-  template <typename FloatType, typename Device>
-  void PrintRayStatus(Ray<FloatType>& rays, Device);
+  template <typename FloatType>
+  void PrintRayStatus(Ray<FloatType>& rays);
 
 protected:
-  template <typename Device, typename FloatType>
-  void RenderOnDevice(Ray<FloatType>& rays, Device);
-
   // Data set info
   vtkm::cont::Field ScalarField;
   vtkm::cont::Field EmissionField;
@@ -180,7 +194,6 @@ protected:
   vtkm::Float32 BoundingBox[6];
 
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>> ColorMap;
-  vtkm::cont::ArrayHandle<vtkm::Id> PreviousCellIds;
 
   vtkm::Vec<vtkm::Float32, 4> BackgroundColor;
   vtkm::Float32 SampleDistance;
@@ -204,6 +217,7 @@ protected:
   vtkm::Float64 SampleTime;
   vtkm::Float64 LostRayTime;
   vtkm::Float64 MeshEntryTime;
+  vtkm::Float32 UnitScalar;
 
 }; // class ConnectivityTracer<CellType,ConnectivityType>
 }
diff --git a/vtkm/rendering/raytracing/ConnectivityTracer.hxx b/vtkm/rendering/raytracing/ConnectivityTracer.hxx
deleted file mode 100644
index ec1aac957..000000000
--- a/vtkm/rendering/raytracing/ConnectivityTracer.hxx
+++ /dev/null
@@ -1,1276 +0,0 @@
-//============================================================================
-//  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 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
-//  Copyright 2015 UT-Battelle, LLC.
-//  Copyright 2015 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/DeviceAdapterAlgorithm.h>
-#include <vtkm/cont/ErrorBadValue.h>
-#include <vtkm/cont/Timer.h>
-#include <vtkm/cont/TryExecute.h>
-#include <vtkm/cont/internal/DeviceAdapterListHelpers.h>
-
-#include <vtkm/rendering/raytracing/Camera.h>
-#include <vtkm/rendering/raytracing/CellIntersector.h>
-#include <vtkm/rendering/raytracing/CellSampler.h>
-#include <vtkm/rendering/raytracing/CellTables.h>
-#include <vtkm/rendering/raytracing/MeshConnectivityBase.h>
-#include <vtkm/rendering/raytracing/Ray.h>
-#include <vtkm/rendering/raytracing/RayOperations.h>
-#include <vtkm/rendering/raytracing/RayTracingTypeDefs.h>
-#include <vtkm/rendering/raytracing/Worklets.h>
-
-#include <vtkm/worklet/DispatcherMapField.h>
-#include <vtkm/worklet/DispatcherMapTopology.h>
-#include <vtkm/worklet/WorkletMapField.h>
-#include <vtkm/worklet/WorkletMapTopology.h>
-
-namespace vtkm
-{
-namespace rendering
-{
-namespace raytracing
-{
-namespace detail
-{
-template <typename FloatType>
-template <typename Device>
-void RayTracking<FloatType>::Compact(vtkm::cont::ArrayHandle<FloatType>& compactedDistances,
-                                     vtkm::cont::ArrayHandle<UInt8>& masks,
-                                     Device)
-{
-  //
-  // These distances are stored in the rays, and it has
-  // already been compacted.
-  //
-  CurrentDistance = compactedDistances;
-
-  vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
-
-  bool distance1IsEnter = EnterDist == &Distance1;
-
-  vtkm::cont::ArrayHandle<FloatType> compactedDistance1;
-  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(Distance1, masks, compactedDistance1);
-  Distance1 = compactedDistance1;
-
-  vtkm::cont::ArrayHandle<FloatType> compactedDistance2;
-  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(Distance2, masks, compactedDistance2);
-  Distance2 = compactedDistance2;
-
-  vtkm::cont::ArrayHandle<vtkm::Int32> compactedExitFace;
-  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(ExitFace, masks, compactedExitFace);
-  ExitFace = compactedExitFace;
-
-  if (distance1IsEnter)
-  {
-    EnterDist = &Distance1;
-    ExitDist = &Distance2;
-  }
-  else
-  {
-    EnterDist = &Distance2;
-    ExitDist = &Distance1;
-  }
-}
-
-template <typename FloatType>
-template <typename Device>
-void RayTracking<FloatType>::Init(const vtkm::Id size,
-                                  vtkm::cont::ArrayHandle<FloatType>& distances,
-                                  Device)
-{
-
-  ExitFace.PrepareForOutput(size, Device());
-  Distance1.PrepareForOutput(size, Device());
-  Distance2.PrepareForOutput(size, Device());
-
-  CurrentDistance = distances;
-  //
-  // Set the initial Distances
-  //
-  vtkm::worklet::DispatcherMapField<CopyAndOffset<FloatType>> resetDistancesDispatcher(
-    CopyAndOffset<FloatType>(0.0f));
-  resetDistancesDispatcher.SetDevice(Device());
-  resetDistancesDispatcher.Invoke(distances, *EnterDist);
-
-  //
-  // Init the exit faces. This value is used to load the next cell
-  // base on the cell and face it left
-  //
-  vtkm::cont::ArrayHandleConstant<vtkm::Int32> negOne(-1, size);
-  vtkm::cont::Algorithm::Copy(Device(), negOne, ExitFace);
-
-  vtkm::cont::ArrayHandleConstant<FloatType> negOnef(-1.f, size);
-  vtkm::cont::Algorithm::Copy(Device(), negOnef, *ExitDist);
-}
-
-template <typename FloatType>
-void RayTracking<FloatType>::Swap()
-{
-  vtkm::cont::ArrayHandle<FloatType>* tmpPtr;
-  tmpPtr = EnterDist;
-  EnterDist = ExitDist;
-  ExitDist = tmpPtr;
-}
-} //namespace detail
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::PrintRayStatus(Ray<FloatType>& rays, Device)
-{
-  vtkm::Id raysExited = RayOperations::GetStatusCount(rays, RAY_EXITED_MESH, Device());
-  vtkm::Id raysActive = RayOperations::GetStatusCount(rays, RAY_ACTIVE, Device());
-  vtkm::Id raysAbandoned = RayOperations::GetStatusCount(rays, RAY_ABANDONED, Device());
-  vtkm::Id raysExitedDom = RayOperations::GetStatusCount(rays, RAY_EXITED_DOMAIN, Device());
-  std::cout << "\r Ray Status " << std::setw(10) << std::left << " Lost " << std::setw(10)
-            << std::left << RaysLost << std::setw(10) << std::left << " Exited " << std::setw(10)
-            << std::left << raysExited << std::setw(10) << std::left << " Active " << std::setw(10)
-            << raysActive << std::setw(10) << std::left << " Abandoned " << std::setw(10)
-            << raysAbandoned << " Exited Domain " << std::setw(10) << std::left << raysExitedDom
-            << "\n";
-}
-
-//
-//  Advance Ray
-//      After a ray leaves the mesh, we need to check to see
-//      of the ray re-enters the mesh within this domain. This
-//      function moves the ray forward some offset to prevent
-//      "shadowing" and hitting the same exit point.
-//
-template <typename FloatType>
-class AdvanceRay : public vtkm::worklet::WorkletMapField
-{
-  FloatType Offset;
-
-public:
-  VTKM_CONT
-  AdvanceRay(const FloatType offset = 0.00001)
-    : Offset(offset)
-  {
-  }
-  using ControlSignature = void(FieldIn<>, FieldInOut<>);
-  using ExecutionSignature = void(_1, _2);
-
-  VTKM_EXEC inline void operator()(const vtkm::UInt8& status, FloatType& distance) const
-  {
-    if (status == RAY_EXITED_MESH)
-      distance += Offset;
-  }
-}; //class AdvanceRay
-
-template <typename FloatType>
-class LocateCell : public vtkm::worklet::WorkletMapField
-{
-private:
-  const MeshConnectivityBase* MeshConn;
-  CellIntersector<255> Intersector;
-
-public:
-  LocateCell(const MeshConnectivityBase* meshConn)
-    : MeshConn(meshConn)
-  {
-  }
-
-  using ControlSignature = void(FieldInOut<>,
-                                WholeArrayIn<>,
-                                FieldIn<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldIn<>);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8);
-
-  template <typename PointPortalType>
-  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
-                                   PointPortalType& vertices,
-                                   const vtkm::Vec<FloatType, 3>& dir,
-                                   FloatType& enterDistance,
-                                   FloatType& exitDistance,
-                                   vtkm::Int32& enterFace,
-                                   vtkm::UInt8& rayStatus,
-                                   const vtkm::Vec<FloatType, 3>& origin) const
-  {
-    if (enterFace != -1 && rayStatus == RAY_ACTIVE)
-    {
-      currentCell = MeshConn->GetConnectingCell(currentCell, enterFace);
-      if (currentCell == -1)
-        rayStatus = RAY_EXITED_MESH;
-      enterFace = -1;
-    }
-    //This ray is dead or exited the mesh and needs re-entry
-    if (rayStatus != RAY_ACTIVE)
-    {
-      return;
-    }
-    FloatType xpoints[8];
-    FloatType ypoints[8];
-    FloatType zpoints[8];
-    vtkm::Id cellConn[8];
-    FloatType distances[6];
-
-    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
-    //load local cell data
-    for (int i = 0; i < numIndices; ++i)
-    {
-      BOUNDS_CHECK(vertices, cellConn[i]);
-      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
-      xpoints[i] = point[0];
-      ypoints[i] = point[1];
-      zpoints[i] = point[2];
-    }
-    const vtkm::UInt8 cellShape = MeshConn->GetCellShape(currentCell);
-    Intersector.IntersectCell(xpoints, ypoints, zpoints, dir, origin, distances, cellShape);
-
-    CellTables tables;
-    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);
-    //vtkm::Int32 minFace = 6;
-    vtkm::Int32 maxFace = -1;
-
-    FloatType minDistance = static_cast<FloatType>(1e32);
-    FloatType maxDistance = static_cast<FloatType>(-1);
-    int hitCount = 0;
-    for (vtkm::Int32 i = 0; i < numFaces; ++i)
-    {
-      FloatType dist = distances[i];
-
-      if (dist != -1)
-      {
-        hitCount++;
-        if (dist < minDistance)
-        {
-          minDistance = dist;
-          //minFace = i;
-        }
-        if (dist > maxDistance)
-        {
-          maxDistance = dist;
-          maxFace = i;
-        }
-      }
-    }
-
-    if (maxDistance <= enterDistance || minDistance == maxDistance)
-    {
-      rayStatus = RAY_LOST;
-    }
-    else
-    {
-      enterDistance = minDistance;
-      exitDistance = maxDistance;
-      enterFace = maxFace;
-    }
-
-  } //operator
-};  //class LocateCell
-
-template <typename FloatType, typename Device>
-class RayBumper : public vtkm::worklet::WorkletMapField
-{
-private:
-  using FloatTypeHandle = typename vtkm::cont::ArrayHandle<FloatType>;
-  using FloatTypePortal = typename FloatTypeHandle::template ExecutionTypes<Device>::Portal;
-  FloatTypePortal DirectionsX;
-  FloatTypePortal DirectionsY;
-  FloatTypePortal DirectionsZ;
-  const MeshConnectivityBase* MeshConn;
-  CellIntersector<255> Intersector;
-  const vtkm::UInt8 FailureStatus; // the status to assign ray if we fail to find the intersection
-public:
-  RayBumper(FloatTypeHandle dirsx,
-            FloatTypeHandle dirsy,
-            FloatTypeHandle dirsz,
-            const MeshConnectivityBase* meshConn,
-            vtkm::UInt8 failureStatus = RAY_ABANDONED)
-    : DirectionsX(dirsx.PrepareForInPlace(Device()))
-    , DirectionsY(dirsy.PrepareForInPlace(Device()))
-    , DirectionsZ(dirsz.PrepareForInPlace(Device()))
-    , MeshConn(meshConn)
-    , FailureStatus(failureStatus)
-  {
-  }
-
-
-  using ControlSignature = void(FieldInOut<>,
-                                WholeArrayIn<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldInOut<>,
-                                FieldIn<>);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex, _6, _7);
-
-  template <typename PointPortalType>
-  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
-                                   PointPortalType& vertices,
-                                   FloatType& enterDistance,
-                                   FloatType& exitDistance,
-                                   vtkm::Int32& enterFace,
-                                   const vtkm::Id& pixelIndex,
-                                   vtkm::UInt8& rayStatus,
-                                   const vtkm::Vec<FloatType, 3>& origin) const
-  {
-    // We only process lost rays
-    if (rayStatus != RAY_LOST)
-    {
-      return;
-    }
-
-    FloatType xpoints[8];
-    FloatType ypoints[8];
-    FloatType zpoints[8];
-    vtkm::Id cellConn[8];
-    FloatType distances[6];
-
-    vtkm::Vec<FloatType, 3> centroid(0., 0., 0.);
-
-    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
-    //load local cell data
-    for (int i = 0; i < numIndices; ++i)
-    {
-      BOUNDS_CHECK(vertices, cellConn[i]);
-      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
-      centroid = centroid + point;
-      xpoints[i] = point[0];
-      ypoints[i] = point[1];
-      zpoints[i] = point[2];
-    }
-
-    FloatType invNumIndices = static_cast<FloatType>(1.) / static_cast<FloatType>(numIndices);
-    centroid[0] = centroid[0] * invNumIndices;
-    centroid[1] = centroid[1] * invNumIndices;
-    centroid[2] = centroid[2] * invNumIndices;
-
-    vtkm::Vec<FloatType, 3> toCentroid = centroid - origin;
-    vtkm::Normalize(toCentroid);
-
-    vtkm::Vec<FloatType, 3> dir(
-      DirectionsX.Get(pixelIndex), DirectionsY.Get(pixelIndex), DirectionsZ.Get(pixelIndex));
-    vtkm::Vec<FloatType, 3> bump = toCentroid - dir;
-    dir = dir + RAY_TUG_EPSILON * bump;
-
-    vtkm::Normalize(dir);
-
-    DirectionsX.Set(pixelIndex, dir[0]);
-    DirectionsY.Set(pixelIndex, dir[1]);
-    DirectionsZ.Set(pixelIndex, dir[2]);
-
-    const vtkm::UInt8 cellShape = MeshConn->GetCellShape(currentCell);
-    Intersector.IntersectCell(xpoints, ypoints, zpoints, dir, origin, distances, cellShape);
-
-    CellTables tables;
-    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);
-
-    //vtkm::Int32 minFace = 6;
-    vtkm::Int32 maxFace = -1;
-    FloatType minDistance = static_cast<FloatType>(1e32);
-    FloatType maxDistance = static_cast<FloatType>(-1);
-    int hitCount = 0;
-    for (int i = 0; i < numFaces; ++i)
-    {
-      FloatType dist = distances[i];
-
-      if (dist != -1)
-      {
-        hitCount++;
-        if (dist < minDistance)
-        {
-          minDistance = dist;
-          //minFace = i;
-        }
-        if (dist >= maxDistance)
-        {
-          maxDistance = dist;
-          maxFace = i;
-        }
-      }
-    }
-    if (minDistance >= maxDistance)
-    {
-      rayStatus = FailureStatus;
-    }
-    else
-    {
-      enterDistance = minDistance;
-      exitDistance = maxDistance;
-      enterFace = maxFace;
-      rayStatus = RAY_ACTIVE; //re-activate ray
-    }
-
-  } //operator
-};  //class RayBumper
-
-template <typename FloatType>
-class AddPathLengths : public vtkm::worklet::WorkletMapField
-{
-public:
-  VTKM_CONT
-  AddPathLengths() {}
-
-  using ControlSignature = void(FieldIn<RayStatusType>,            // ray status
-                                FieldIn<ScalarRenderingTypes>,     // cell enter distance
-                                FieldIn<ScalarRenderingTypes>,     // cell exit distance
-                                FieldInOut<ScalarRenderingTypes>); // ray absorption data
-
-  using ExecutionSignature = void(_1, _2, _3, _4);
-
-  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
-                                   const FloatType& enterDistance,
-                                   const FloatType& exitDistance,
-                                   FloatType& distance) const
-  {
-    if (rayStatus != RAY_ACTIVE)
-    {
-      return;
-    }
-
-    if (exitDistance <= enterDistance)
-    {
-      return;
-    }
-
-    FloatType segmentLength = exitDistance - enterDistance;
-    distance += segmentLength;
-  }
-};
-
-template <typename FloatType>
-class Integrate : public vtkm::worklet::WorkletMapField
-{
-private:
-  const vtkm::Int32 NumBins;
-
-public:
-  VTKM_CONT
-  Integrate(const vtkm::Int32 numBins)
-    : NumBins(numBins)
-  {
-  }
-
-  using ControlSignature = void(FieldIn<RayStatusType>,                // ray status
-                                FieldIn<ScalarRenderingTypes>,         // cell enter distance
-                                FieldIn<ScalarRenderingTypes>,         // cell exit distance
-                                FieldInOut<ScalarRenderingTypes>,      // current distance
-                                WholeArrayIn<ScalarRenderingTypes>,    // cell absorption data array
-                                WholeArrayInOut<ScalarRenderingTypes>, // ray absorption data
-                                FieldIn<IdType>);                      // current cell
-
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, WorkIndex);
-
-  template <typename CellDataPortalType, typename RayDataPortalType>
-  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
-                                   const FloatType& enterDistance,
-                                   const FloatType& exitDistance,
-                                   FloatType& currentDistance,
-                                   const CellDataPortalType& cellData,
-                                   RayDataPortalType& energyBins,
-                                   const vtkm::Id& currentCell,
-                                   const vtkm::Id& rayIndex) const
-  {
-    if (rayStatus != RAY_ACTIVE)
-    {
-      return;
-    }
-    if (exitDistance <= enterDistance)
-    {
-      return;
-    }
-
-    FloatType segmentLength = exitDistance - enterDistance;
-
-    vtkm::Id rayOffset = NumBins * rayIndex;
-    vtkm::Id cellOffset = NumBins * currentCell;
-    for (vtkm::Int32 i = 0; i < NumBins; ++i)
-    {
-      BOUNDS_CHECK(cellData, cellOffset + i);
-      FloatType absorb = static_cast<FloatType>(cellData.Get(cellOffset + i));
-
-      absorb = vtkm::Exp(-absorb * segmentLength);
-      BOUNDS_CHECK(energyBins, rayOffset + i);
-      FloatType intensity = static_cast<FloatType>(energyBins.Get(rayOffset + i));
-      energyBins.Set(rayOffset + i, intensity * absorb);
-    }
-    currentDistance = exitDistance;
-  }
-};
-
-template <typename FloatType>
-class IntegrateEmission : public vtkm::worklet::WorkletMapField
-{
-private:
-  const vtkm::Int32 NumBins;
-  bool DivideEmisByAbsorb;
-
-public:
-  VTKM_CONT
-  IntegrateEmission(const vtkm::Int32 numBins, const bool divideEmisByAbsorb)
-    : NumBins(numBins)
-    , DivideEmisByAbsorb(divideEmisByAbsorb)
-  {
-  }
-
-  using ControlSignature = void(FieldIn<>,                          // ray status
-                                FieldIn<>,                          // cell enter distance
-                                FieldIn<>,                          // cell exit distance
-                                FieldInOut<>,                       // current distance
-                                WholeArrayIn<ScalarRenderingTypes>, // cell absorption data array
-                                WholeArrayIn<ScalarRenderingTypes>, // cell emission data array
-                                WholeArrayInOut<>,                  // ray absorption data
-                                WholeArrayInOut<>,                  // ray emission data
-                                FieldIn<>);                         // current cell
-
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9, WorkIndex);
-
-  template <typename CellAbsPortalType, typename CellEmisPortalType, typename RayDataPortalType>
-  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
-                                   const FloatType& enterDistance,
-                                   const FloatType& exitDistance,
-                                   FloatType& currentDistance,
-                                   const CellAbsPortalType& absorptionData,
-                                   const CellEmisPortalType& emissionData,
-                                   RayDataPortalType& absorptionBins,
-                                   RayDataPortalType& emissionBins,
-                                   const vtkm::Id& currentCell,
-                                   const vtkm::Id& rayIndex) const
-  {
-    if (rayStatus != RAY_ACTIVE)
-    {
-      return;
-    }
-    if (exitDistance <= enterDistance)
-    {
-      return;
-    }
-
-    FloatType segmentLength = exitDistance - enterDistance;
-
-    vtkm::Id rayOffset = NumBins * rayIndex;
-    vtkm::Id cellOffset = NumBins * currentCell;
-    for (vtkm::Int32 i = 0; i < NumBins; ++i)
-    {
-      BOUNDS_CHECK(absorptionData, cellOffset + i);
-      FloatType absorb = static_cast<FloatType>(absorptionData.Get(cellOffset + i));
-      BOUNDS_CHECK(emissionData, cellOffset + i);
-      FloatType emission = static_cast<FloatType>(emissionData.Get(cellOffset + i));
-
-      if (DivideEmisByAbsorb)
-      {
-        emission /= absorb;
-      }
-
-      FloatType tmp = vtkm::Exp(-absorb * segmentLength);
-      BOUNDS_CHECK(absorptionBins, rayOffset + i);
-
-      //
-      // Traditionally, we would only keep track of a single intensity value per ray
-      // per bin and we would integrate from the beginning to end of the ray. In a
-      // distributed memory setting, we would move cell data around so that the
-      // entire ray could be traced, but in situ, moving that much cell data around
-      // could blow memory. Here we are keeping track of two values. Total absorption
-      // through this contiguous segment of the mesh, and the amount of emitted energy
-      // that makes it out of this mesh segment. If this is really run on a single node,
-      // we can get the final energy value by multiplying the background intensity by
-      // the total absorption of the mesh segment and add in the amount of emitted
-      // energy that escapes.
-      //
-      FloatType absorbIntensity = static_cast<FloatType>(absorptionBins.Get(rayOffset + i));
-      FloatType emissionIntensity = static_cast<FloatType>(emissionBins.Get(rayOffset + i));
-
-      absorptionBins.Set(rayOffset + i, absorbIntensity * tmp);
-
-      emissionIntensity = emissionIntensity * tmp + emission * (1.f - tmp);
-
-      BOUNDS_CHECK(emissionBins, rayOffset + i);
-      emissionBins.Set(rayOffset + i, emissionIntensity);
-    }
-    currentDistance = exitDistance;
-  }
-};
-//
-//  IdentifyMissedRay is a debugging routine that detects
-//  rays that fail to have any value because of a external
-//  intersection and cell intersection mismatch
-//
-//
-class IdentifyMissedRay : public vtkm::worklet::WorkletMapField
-{
-public:
-  vtkm::Id Width;
-  vtkm::Id Height;
-  vtkm::Vec<vtkm::Float32, 4> BGColor;
-  IdentifyMissedRay(const vtkm::Id width,
-                    const vtkm::Id height,
-                    vtkm::Vec<vtkm::Float32, 4> bgcolor)
-    : Width(width)
-    , Height(height)
-    , BGColor(bgcolor)
-  {
-  }
-  using ControlSignature = void(FieldIn<>, WholeArrayIn<>);
-  using ExecutionSignature = void(_1, _2);
-
-
-  VTKM_EXEC inline bool IsBGColor(const vtkm::Vec<vtkm::Float32, 4> color) const
-  {
-    bool isBG = false;
-
-    if (color[0] == BGColor[0] && color[1] == BGColor[1] && color[2] == BGColor[2] &&
-        color[3] == BGColor[3])
-      isBG = true;
-    return isBG;
-  }
-
-  template <typename ColorBufferType>
-  VTKM_EXEC inline void operator()(const vtkm::Id& pixelId, ColorBufferType& buffer) const
-  {
-    vtkm::Id x = pixelId % Width;
-    vtkm::Id y = pixelId / Width;
-
-    // Conservative check, we only want to check pixels in the middle
-    if (x <= 0 || y <= 0)
-      return;
-    if (x >= Width - 1 || y >= Height - 1)
-      return;
-    vtkm::Vec<vtkm::Float32, 4> pixel;
-    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 0));
-    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 1));
-    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 2));
-    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 3));
-    if (!IsBGColor(pixel))
-      return;
-    vtkm::Id p0 = (y)*Width + (x + 1);
-    vtkm::Id p1 = (y)*Width + (x - 1);
-    vtkm::Id p2 = (y + 1) * Width + (x);
-    vtkm::Id p3 = (y - 1) * Width + (x);
-    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 0));
-    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 1));
-    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 2));
-    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 3));
-    if (IsBGColor(pixel))
-      return;
-    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 0));
-    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 1));
-    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 2));
-    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 3));
-    if (IsBGColor(pixel))
-      return;
-    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 0));
-    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 1));
-    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 2));
-    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 3));
-    if (IsBGColor(pixel))
-      return;
-    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 0));
-    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 1));
-    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 2));
-    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 3));
-    if (IsBGColor(pixel))
-      return;
-
-    printf("Possible error ray missed ray %d\n", (int)pixelId);
-  }
-};
-
-template <typename FloatType, typename Device>
-class SampleCellAssocCells : public vtkm::worklet::WorkletMapField
-{
-private:
-  using ColorHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
-  using ColorBuffer = typename vtkm::cont::ArrayHandle<FloatType>;
-  using ColorConstPortal = typename ColorHandle::ExecutionTypes<Device>::PortalConst;
-  using ColorPortal = typename ColorBuffer::template ExecutionTypes<Device>::Portal;
-
-  CellSampler<255> Sampler;
-  FloatType SampleDistance;
-  FloatType MinScalar;
-  FloatType InvDeltaScalar;
-  ColorPortal FrameBuffer;
-  ColorConstPortal ColorMap;
-  const MeshConnectivityBase* MeshConn;
-  vtkm::Int32 ColorMapSize;
-
-public:
-  SampleCellAssocCells(const FloatType& sampleDistance,
-                       const FloatType& minScalar,
-                       const FloatType& maxScalar,
-                       ColorHandle& colorMap,
-                       ColorBuffer& frameBuffer,
-                       const MeshConnectivityBase* meshConn)
-    : SampleDistance(sampleDistance)
-    , MinScalar(minScalar)
-    , ColorMap(colorMap.PrepareForInput(Device()))
-    , MeshConn(meshConn)
-  {
-    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
-    ColorMapSize = static_cast<vtkm::Int32>(ColorMap.GetNumberOfValues());
-    this->FrameBuffer = frameBuffer.PrepareForOutput(frameBuffer.GetNumberOfValues(), Device());
-  }
-
-
-  using ControlSignature = void(FieldIn<>,
-                                WholeArrayIn<ScalarRenderingTypes>,
-                                FieldIn<>,
-                                FieldIn<>,
-                                FieldInOut<>,
-                                FieldInOut<>);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex);
-
-  template <typename ScalarPortalType>
-  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
-                                   ScalarPortalType& scalarPortal,
-                                   const FloatType& enterDistance,
-                                   const FloatType& exitDistance,
-                                   FloatType& currentDistance,
-                                   vtkm::UInt8& rayStatus,
-                                   const vtkm::Id& pixelIndex) const
-  {
-
-    if (rayStatus != RAY_ACTIVE)
-      return;
-
-    vtkm::Vec<vtkm::Float32, 4> color;
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
-    color[0] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 0));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
-    color[1] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 1));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
-    color[2] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 2));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
-    color[3] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 3));
-
-    vtkm::Float32 scalar;
-    BOUNDS_CHECK(scalarPortal, currentCell);
-    scalar = vtkm::Float32(scalarPortal.Get(currentCell));
-    //
-    // There can be mismatches in the initial enter distance and the current distance
-    // due to lost rays at cell borders. For now,
-    // we will just advance the current position to the enter distance, since otherwise,
-    // the pixel would never be sampled.
-    //
-    if (currentDistance < enterDistance)
-      currentDistance = enterDistance;
-
-    vtkm::Float32 lerpedScalar;
-    lerpedScalar = static_cast<vtkm::Float32>((scalar - MinScalar) * InvDeltaScalar);
-    vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(ColorMapSize));
-    if (colorIndex < 0)
-      colorIndex = 0;
-    if (colorIndex >= ColorMapSize)
-      colorIndex = ColorMapSize - 1;
-    BOUNDS_CHECK(ColorMap, colorIndex);
-    vtkm::Vec<vtkm::Float32, 4> sampleColor = ColorMap.Get(colorIndex);
-
-    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
-    {
-      //composite
-      vtkm::Float32 alpha = sampleColor[3] * (1.f - color[3]);
-      color[0] = color[0] + sampleColor[0] * alpha;
-      color[1] = color[1] + sampleColor[1] * alpha;
-      color[2] = color[2] + sampleColor[2] * alpha;
-      color[3] = alpha + color[3];
-
-      if (color[3] > 1.)
-      {
-        rayStatus = RAY_TERMINATED;
-        break;
-      }
-      currentDistance += SampleDistance;
-    }
-
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
-    FrameBuffer.Set(pixelIndex * 4 + 0, color[0]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
-    FrameBuffer.Set(pixelIndex * 4 + 1, color[1]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
-    FrameBuffer.Set(pixelIndex * 4 + 2, color[2]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
-    FrameBuffer.Set(pixelIndex * 4 + 3, color[3]);
-  }
-}; //class Sample cell
-
-template <typename FloatType, typename Device>
-class SampleCellAssocPoints : public vtkm::worklet::WorkletMapField
-{
-private:
-  using ColorHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
-  using ColorBuffer = typename vtkm::cont::ArrayHandle<FloatType>;
-  using ColorConstPortal = typename ColorHandle::ExecutionTypes<Device>::PortalConst;
-  using ColorPortal = typename ColorBuffer::template ExecutionTypes<Device>::Portal;
-
-  CellSampler<255> Sampler;
-  FloatType SampleDistance;
-  const MeshConnectivityBase* MeshConn;
-  FloatType MinScalar;
-  FloatType InvDeltaScalar;
-  ColorPortal FrameBuffer;
-  ColorConstPortal ColorMap;
-  vtkm::Id ColorMapSize;
-
-public:
-  SampleCellAssocPoints(const FloatType& sampleDistance,
-                        const FloatType& minScalar,
-                        const FloatType& maxScalar,
-                        ColorHandle& colorMap,
-                        ColorBuffer& frameBuffer,
-                        const MeshConnectivityBase* meshConn)
-    : SampleDistance(sampleDistance)
-    , MeshConn(meshConn)
-    , MinScalar(minScalar)
-    , ColorMap(colorMap.PrepareForInput(Device()))
-  {
-    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
-    ColorMapSize = ColorMap.GetNumberOfValues();
-    this->FrameBuffer = frameBuffer.PrepareForOutput(frameBuffer.GetNumberOfValues(), Device());
-  }
-
-
-  using ControlSignature = void(FieldIn<>,
-                                WholeArrayIn<Vec3>,
-                                WholeArrayIn<ScalarRenderingTypes>,
-                                FieldIn<>,
-                                FieldIn<>,
-                                FieldInOut<>,
-                                FieldIn<>,
-                                FieldInOut<>,
-                                FieldIn<>);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, WorkIndex, _9);
-
-  template <typename PointPortalType, typename ScalarPortalType>
-  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
-                                   PointPortalType& vertices,
-                                   ScalarPortalType& scalarPortal,
-                                   const FloatType& enterDistance,
-                                   const FloatType& exitDistance,
-                                   FloatType& currentDistance,
-                                   const vtkm::Vec<vtkm::Float32, 3>& dir,
-                                   vtkm::UInt8& rayStatus,
-                                   const vtkm::Id& pixelIndex,
-                                   const vtkm::Vec<FloatType, 3>& origin) const
-  {
-
-    if (rayStatus != RAY_ACTIVE)
-      return;
-
-    vtkm::Vec<vtkm::Float32, 4> color;
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
-    color[0] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 0));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
-    color[1] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 1));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
-    color[2] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 2));
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
-    color[3] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 3));
-
-    if (color[3] >= 1.f)
-    {
-      rayStatus = RAY_TERMINATED;
-      return;
-    }
-    vtkm::Vec<vtkm::Float32, 8> scalars;
-    vtkm::Vec<vtkm::Vec<FloatType, 3>, 8> points;
-    // silence "may" be uninitialized warning
-    for (vtkm::Int32 i = 0; i < 8; ++i)
-    {
-      scalars[i] = 0.f;
-      points[i] = vtkm::Vec<FloatType, 3>(0.f, 0.f, 0.f);
-    }
-    //load local scalar cell data
-    vtkm::Id cellConn[8];
-    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
-    for (int i = 0; i < numIndices; ++i)
-    {
-      BOUNDS_CHECK(scalarPortal, cellConn[i]);
-      scalars[i] = static_cast<vtkm::Float32>(scalarPortal.Get(cellConn[i]));
-      BOUNDS_CHECK(vertices, cellConn[i]);
-      points[i] = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
-    }
-    //
-    // There can be mismatches in the initial enter distance and the current distance
-    // due to lost rays at cell borders. For now,
-    // we will just advance the current position to the enter distance, since otherwise,
-    // the pixel would never be sampled.
-    //
-    if (currentDistance < enterDistance)
-    {
-      currentDistance = enterDistance;
-    }
-
-    const vtkm::Int32 cellShape = MeshConn->GetCellShape(currentCell);
-    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
-    {
-      vtkm::Vec<FloatType, 3> sampleLoc = origin + currentDistance * dir;
-      vtkm::Float32 lerpedScalar;
-      bool validSample =
-        Sampler.SampleCell(points, scalars, sampleLoc, lerpedScalar, *this, cellShape);
-      if (!validSample)
-      {
-        //
-        // There is a slight mismatch between intersections and parametric coordinates
-        // which results in a invalid sample very close to the cell edge. Just throw
-        // this sample away, and move to the next sample.
-        //
-
-        //There should be a sample here, so offset and try again.
-
-        currentDistance += 0.00001f;
-        continue;
-      }
-      lerpedScalar = static_cast<vtkm::Float32>((lerpedScalar - MinScalar) * InvDeltaScalar);
-      vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(ColorMapSize));
-
-      colorIndex = vtkm::Min(vtkm::Max(colorIndex, vtkm::Id(0)), ColorMapSize - 1);
-      BOUNDS_CHECK(ColorMap, colorIndex);
-      vtkm::Vec<vtkm::Float32, 4> sampleColor = ColorMap.Get(colorIndex);
-      //composite
-      sampleColor[3] *= (1.f - color[3]);
-      color[0] = color[0] + sampleColor[0] * sampleColor[3];
-      color[1] = color[1] + sampleColor[1] * sampleColor[3];
-      color[2] = color[2] + sampleColor[2] * sampleColor[3];
-      color[3] = sampleColor[3] + color[3];
-
-      if (color[3] >= 1.0)
-      {
-        rayStatus = RAY_TERMINATED;
-        break;
-      }
-      currentDistance += SampleDistance;
-    }
-
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
-    FrameBuffer.Set(pixelIndex * 4 + 0, color[0]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
-    FrameBuffer.Set(pixelIndex * 4 + 1, color[1]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
-    FrameBuffer.Set(pixelIndex * 4 + 2, color[2]);
-    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
-    FrameBuffer.Set(pixelIndex * 4 + 3, color[3]);
-  }
-}; //class Sample cell
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::IntersectCell(Ray<FloatType>& rays,
-                                       detail::RayTracking<FloatType>& tracker,
-                                       const MeshConnectivityBase* meshConn,
-                                       Device)
-{
-  using LocateC = LocateCell<FloatType>;
-  vtkm::cont::Timer<Device> timer;
-  vtkm::worklet::DispatcherMapField<LocateC> locateDispatch(meshConn);
-  locateDispatch.SetDevice(Device());
-  locateDispatch.Invoke(rays.HitIdx,
-                        this->Coords,
-                        rays.Dir,
-                        *(tracker.EnterDist),
-                        *(tracker.ExitDist),
-                        tracker.ExitFace,
-                        rays.Status,
-                        rays.Origin);
-
-  if (this->CountRayStatus)
-    RaysLost = RayOperations::GetStatusCount(rays, RAY_LOST, Device());
-  this->IntersectTime += timer.GetElapsedTime();
-}
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::AccumulatePathLengths(Ray<FloatType>& rays,
-                                               detail::RayTracking<FloatType>& tracker,
-                                               Device)
-{
-  vtkm::worklet::DispatcherMapField<AddPathLengths<FloatType>> dispatcher(
-    (AddPathLengths<FloatType>()));
-  dispatcher.SetDevice(Device());
-  dispatcher.Invoke(
-    rays.Status, *(tracker.EnterDist), *(tracker.ExitDist), rays.GetBuffer("path_lengths").Buffer);
-}
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::FindLostRays(Ray<FloatType>& rays,
-                                      detail::RayTracking<FloatType>& tracker,
-                                      const MeshConnectivityBase* meshConn,
-                                      Device)
-{
-  using RayB = RayBumper<FloatType, Device>;
-  vtkm::cont::Timer<Device> timer;
-
-  vtkm::worklet::DispatcherMapField<RayB> bumpDispatch(
-    RayB(rays.DirX, rays.DirY, rays.DirZ, meshConn));
-  bumpDispatch.SetDevice(Device());
-  bumpDispatch.Invoke(rays.HitIdx,
-                      this->Coords,
-                      *(tracker.EnterDist),
-                      *(tracker.ExitDist),
-                      tracker.ExitFace,
-                      rays.Status,
-                      rays.Origin);
-
-  this->LostRayTime += timer.GetElapsedTime();
-}
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::SampleCells(Ray<FloatType>& rays,
-                                     detail::RayTracking<FloatType>& tracker,
-                                     const MeshConnectivityBase* meshConn,
-                                     Device)
-{
-  using SampleP = SampleCellAssocPoints<FloatType, Device>;
-  using SampleC = SampleCellAssocCells<FloatType, Device>;
-  vtkm::cont::Timer<Device> timer;
-
-  VTKM_ASSERT(rays.Buffers.at(0).GetNumChannels() == 4);
-
-  if (FieldAssocPoints)
-  {
-    vtkm::worklet::DispatcherMapField<SampleP> dispatcher(
-      SampleP(this->SampleDistance,
-              vtkm::Float32(this->ScalarBounds.Min),
-              vtkm::Float32(this->ScalarBounds.Max),
-              this->ColorMap,
-              rays.Buffers.at(0).Buffer,
-              meshConn));
-    dispatcher.SetDevice(Device());
-    dispatcher.Invoke(rays.HitIdx,
-                      this->Coords,
-                      this->ScalarField.GetData(),
-                      *(tracker.EnterDist),
-                      *(tracker.ExitDist),
-                      tracker.CurrentDistance,
-                      rays.Dir,
-                      rays.Status,
-                      rays.Origin);
-  }
-  else
-  {
-    vtkm::worklet::DispatcherMapField<SampleC> dispatcher(
-      SampleC(this->SampleDistance,
-              vtkm::Float32(this->ScalarBounds.Min),
-              vtkm::Float32(this->ScalarBounds.Max),
-              this->ColorMap,
-              rays.Buffers.at(0).Buffer,
-              meshConn));
-    dispatcher.SetDevice(Device());
-    dispatcher.Invoke(rays.HitIdx,
-                      this->ScalarField.GetData(),
-                      *(tracker.EnterDist),
-                      *(tracker.ExitDist),
-                      tracker.CurrentDistance,
-                      rays.Status);
-  }
-
-  this->SampleTime += timer.GetElapsedTime();
-}
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::IntegrateCells(Ray<FloatType>& rays,
-                                        detail::RayTracking<FloatType>& tracker,
-                                        Device)
-{
-  vtkm::cont::Timer<Device> timer;
-  if (HasEmission)
-  {
-    bool divideEmisByAbsorp = false;
-    vtkm::cont::ArrayHandle<FloatType> absorp = rays.Buffers.at(0).Buffer;
-    vtkm::cont::ArrayHandle<FloatType> emission = rays.GetBuffer("emission").Buffer;
-    vtkm::worklet::DispatcherMapField<IntegrateEmission<FloatType>> dispatcher(
-      IntegrateEmission<FloatType>(rays.Buffers.at(0).GetNumChannels(), divideEmisByAbsorp));
-    dispatcher.SetDevice(Device());
-    dispatcher.Invoke(rays.Status,
-                      *(tracker.EnterDist),
-                      *(tracker.ExitDist),
-                      rays.Distance,
-                      this->ScalarField.GetData(),
-                      this->EmissionField.GetData(),
-                      absorp,
-                      emission,
-                      rays.HitIdx);
-  }
-  else
-  {
-    vtkm::worklet::DispatcherMapField<Integrate<FloatType>> dispatcher(
-      Integrate<FloatType>(rays.Buffers.at(0).GetNumChannels()));
-    dispatcher.SetDevice(Device());
-    dispatcher.Invoke(rays.Status,
-                      *(tracker.EnterDist),
-                      *(tracker.ExitDist),
-                      rays.Distance,
-                      this->ScalarField.GetData(),
-                      rays.Buffers.at(0).Buffer,
-                      rays.HitIdx);
-  }
-
-  IntegrateTime += timer.GetElapsedTime();
-}
-
-// template <vtkm::Int32 CellType>
-// template <typename FloatType>
-// void ConnectivityTracer<CellType>::PrintDebugRay(Ray<FloatType>& rays, vtkm::Id rayId)
-// {
-//   vtkm::Id index = -1;
-//   for (vtkm::Id i = 0; i < rays.NumRays; ++i)
-//   {
-//     if (rays.PixelIdx.GetPortalControl().Get(i) == rayId)
-//     {
-//       index = i;
-//       break;
-//     }
-//   }
-//   if (index == -1)
-//   {
-//     return;
-//   }
-
-//   std::cout << "++++++++RAY " << rayId << "++++++++\n";
-//   std::cout << "Status: " << (int)rays.Status.GetPortalControl().Get(index) << "\n";
-//   std::cout << "HitIndex: " << rays.HitIdx.GetPortalControl().Get(index) << "\n";
-//   std::cout << "Dist " << rays.Distance.GetPortalControl().Get(index) << "\n";
-//   std::cout << "MinDist " << rays.MinDistance.GetPortalControl().Get(index) << "\n";
-//   std::cout << "Origin " << rays.Origin.GetPortalConstControl().Get(index) << "\n";
-//   std::cout << "Dir " << rays.Dir.GetPortalConstControl().Get(index) << "\n";
-//   std::cout << "+++++++++++++++++++++++++\n";
-// }
-
-template <typename FloatType, typename Device>
-void ConnectivityTracer::OffsetMinDistances(Ray<FloatType>& rays, Device)
-{
-  vtkm::worklet::DispatcherMapField<AdvanceRay<FloatType>> dispatcher(
-    AdvanceRay<FloatType>(FloatType(0.001)));
-  dispatcher.SetDevice(Device());
-  dispatcher.Invoke(rays.Status, rays.MinDistance);
-}
-
-template <typename Device, typename FloatType>
-void ConnectivityTracer::RenderOnDevice(Ray<FloatType>& rays, Device)
-{
-  this->RaysLost = 0;
-  Logger* logger = Logger::GetInstance();
-  logger->OpenLogEntry("conn_tracer");
-  logger->AddLogData("device", GetDeviceString(Device()));
-  this->ResetTimers();
-  vtkm::cont::Timer<Device> renderTimer;
-
-  vtkm::cont::DeviceAdapterId devId = Device();
-
-  const MeshConnectivityBase* meshConn = MeshContainer->Construct(devId);
-
-  bool hasPathLengths = rays.HasBuffer("path_lengths");
-
-  vtkm::cont::Timer<Device> timer;
-  this->Init(); // sets sample distance
-  //
-  // All Rays begin as exited to force intersection
-  //
-  RayOperations::ResetStatus(rays, RAY_EXITED_MESH, Device());
-
-  detail::RayTracking<FloatType> rayTracker;
-
-  rayTracker.Init(rays.NumRays, rays.Distance, Device());
-  vtkm::Float64 time = timer.GetElapsedTime();
-  logger->AddLogData("init", time);
-
-  //CountRayStatus = true;
-  bool cullMissedRays = true;
-  bool workRemaining = true;
-  if (this->CountRayStatus)
-  {
-    this->PrintRayStatus(rays, Device());
-  }
-
-  do
-  {
-    {
-      vtkm::cont::Timer<Device> entryTimer;
-      //
-      // if ray misses the external face it will be marked RAY_EXITED_MESH
-      //
-      MeshContainer->FindEntry(rays, devId);
-      MeshEntryTime += entryTimer.GetElapsedTime();
-    }
-
-    if (this->CountRayStatus)
-    {
-      this->PrintRayStatus(rays, Device());
-    }
-    if (cullMissedRays)
-    {
-      //TODO: if we always call this after intersection, then
-      //      we could make a specialized version that only compacts
-      //      hitIdx distance and status, resizing everything else.
-      vtkm::cont::ArrayHandle<UInt8> activeRays;
-      activeRays = RayOperations::CompactActiveRays(rays, Device());
-      rayTracker.Compact(rays.Distance, activeRays, Device());
-      cullMissedRays = false;
-    }
-
-    if (this->CountRayStatus)
-    {
-      this->PrintRayStatus(rays, Device());
-    }
-    // TODO: we should compact out exited rays once below a threshold
-    while (RayOperations::RaysInMesh(rays, Device()))
-    {
-      //
-      // Rays the leave the mesh will be marked as RAYEXITED_MESH
-      this->IntersectCell(rays, rayTracker, meshConn, Device());
-      //
-      // If the ray was lost due to precision issues, we find it.
-      // If it is marked RAY_ABANDONED, then something went wrong.
-      //
-      this->FindLostRays(rays, rayTracker, meshConn, Device());
-      //
-      // integrate along the ray
-      //
-      if (this->Integrator == Volume)
-        this->SampleCells(rays, rayTracker, meshConn, Device());
-      else
-        this->IntegrateCells(rays, rayTracker, Device());
-
-      if (hasPathLengths)
-      {
-        this->AccumulatePathLengths(rays, rayTracker, Device());
-      }
-      //swap enter and exit distances
-      rayTracker.Swap();
-      if (this->CountRayStatus)
-        this->PrintRayStatus(rays, Device());
-    } //for
-
-    workRemaining = RayOperations::RaysProcessed(rays, Device()) != rays.NumRays;
-    //
-    // Ensure that we move the current distance forward some
-    // epsilon so we don't re-enter the cell we just left.
-    //
-    if (workRemaining)
-    {
-      RayOperations::CopyDistancesToMin(rays, Device());
-      this->OffsetMinDistances(rays, Device());
-    }
-  } while (workRemaining);
-
-  if (rays.DebugWidth != -1 && this->Integrator == Volume)
-  {
-
-    vtkm::cont::ArrayHandleCounting<vtkm::Id> pCounter(0, 1, rays.NumRays);
-    vtkm::worklet::DispatcherMapField<IdentifyMissedRay> dispatcher(
-      IdentifyMissedRay(rays.DebugWidth, rays.DebugHeight, this->BackgroundColor));
-    dispatcher.SetDevice(Device());
-    dispatcher.Invoke(pCounter, rays.Buffers.at(0).Buffer);
-  }
-  vtkm::Float64 renderTime = renderTimer.GetElapsedTime();
-  this->LogTimers();
-  logger->AddLogData("active_pixels", rays.NumRays);
-  logger->CloseLogEntry(renderTime);
-} //Render
-}
-}
-} // namespace vtkm::rendering::raytracing
diff --git a/vtkm/rendering/raytracing/CylinderIntersector.cxx b/vtkm/rendering/raytracing/CylinderIntersector.cxx
index 82b5c1545..a3a912fc3 100644
--- a/vtkm/rendering/raytracing/CylinderIntersector.cxx
+++ b/vtkm/rendering/raytracing/CylinderIntersector.cxx
@@ -34,8 +34,126 @@ namespace raytracing
 namespace detail
 {
 
+class FindCylinderAABBs : public vtkm::worklet::WorkletMapField
+{
+public:
+  VTKM_CONT
+  FindCylinderAABBs() {}
+  typedef void ControlSignature(FieldIn<>,
+                                FieldIn<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                WholeArrayIn<Vec3RenderingTypes>);
+  typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8, _9);
+  template <typename PointPortalType>
+  VTKM_EXEC void operator()(const vtkm::Id3 cylId,
+                            const vtkm::Float32& radius,
+                            vtkm::Float32& xmin,
+                            vtkm::Float32& ymin,
+                            vtkm::Float32& zmin,
+                            vtkm::Float32& xmax,
+                            vtkm::Float32& ymax,
+                            vtkm::Float32& zmax,
+                            const PointPortalType& points) const
+  {
+    // cast to Float32
+    vtkm::Vec<vtkm::Float32, 3> point1, point2;
+    vtkm::Vec<vtkm::Float32, 3> temp;
+
+    point1 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[1]));
+    point2 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[2]));
+
+    temp[0] = radius;
+    temp[1] = 0.0f;
+    temp[2] = 0.0f;
+    xmin = ymin = zmin = vtkm::Infinity32();
+    xmax = ymax = zmax = vtkm::NegativeInfinity32();
+
+
+    //set first point to max and min
+    Bounds(point1, radius, xmin, ymin, zmin, xmax, ymax, zmax);
+
+    Bounds(point2, radius, xmin, ymin, zmin, xmax, ymax, zmax);
+  }
+
+  VTKM_EXEC void Bounds(const vtkm::Vec<vtkm::Float32, 3>& point,
+                        const vtkm::Float32& radius,
+                        vtkm::Float32& xmin,
+                        vtkm::Float32& ymin,
+                        vtkm::Float32& zmin,
+                        vtkm::Float32& xmax,
+                        vtkm::Float32& ymax,
+                        vtkm::Float32& zmax) const
+  {
+    vtkm::Vec<vtkm::Float32, 3> temp, p;
+    temp[0] = radius;
+    temp[1] = 0.0f;
+    temp[2] = 0.0f;
+    p = point + temp;
+
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+
+    p = point - temp;
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+
+    temp[0] = 0.f;
+    temp[1] = radius;
+    temp[2] = 0.f;
+
+    p = point + temp;
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+
+    p = point - temp;
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+
+    temp[0] = 0.f;
+    temp[1] = 0.f;
+    temp[2] = radius;
+
+    p = point + temp;
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+
+    p = point - temp;
+    xmin = vtkm::Min(xmin, p[0]);
+    xmax = vtkm::Max(xmax, p[0]);
+    ymin = vtkm::Min(ymin, p[1]);
+    ymax = vtkm::Max(ymax, p[1]);
+    zmin = vtkm::Min(zmin, p[2]);
+    zmax = vtkm::Max(zmax, p[2]);
+  }
+}; //class FindAABBs
+
 template <typename Device>
-class CylinderLeafIntersector : public vtkm::cont::ExecutionObjectBase
+class CylinderLeafIntersector
 {
 public:
   using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id3>;
@@ -45,6 +163,8 @@ public:
   IdArrayPortal CylIds;
   FloatPortal Radii;
 
+  CylinderLeafIntersector() {}
+
   CylinderLeafIntersector(const IdHandle& cylIds, const FloatHandle& radii)
     : CylIds(cylIds.PrepareForInput(Device()))
     , Radii(radii.PrepareForInput(Device()))
@@ -183,17 +303,25 @@ public:
   }
 };
 
-struct IntersectFunctor
+class CylinderLeafWrapper : public vtkm::cont::ExecutionObjectBase
 {
-  template <typename Device, typename Precision>
-  VTKM_CONT bool operator()(Device,
-                            CylinderIntersector* self,
-                            Ray<Precision>& rays,
-                            bool returnCellIndex)
+protected:
+  using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id3>;
+  using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
+  IdHandle CylIds;
+  FloatHandle Radii;
+
+public:
+  CylinderLeafWrapper(IdHandle& cylIds, FloatHandle radii)
+    : CylIds(cylIds)
+    , Radii(radii)
   {
-    VTKM_IS_DEVICE_ADAPTER_TAG(Device);
-    self->IntersectRaysImp(Device(), rays, returnCellIndex);
-    return true;
+  }
+
+  template <typename Device>
+  VTKM_CONT CylinderLeafIntersector<Device> PrepareForExecution(Device) const
+  {
+    return CylinderLeafIntersector<Device>(CylIds, Radii);
   }
 };
 
@@ -302,27 +430,47 @@ CylinderIntersector::~CylinderIntersector()
 {
 }
 
+void CylinderIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
+                                  vtkm::cont::ArrayHandle<vtkm::Id3> cylIds,
+                                  vtkm::cont::ArrayHandle<vtkm::Float32> radii)
+{
+  this->Radii = radii;
+  this->CylIds = cylIds;
+  this->CoordsHandle = coords;
+  AABBs AABB;
+
+  vtkm::worklet::DispatcherMapField<detail::FindCylinderAABBs>(detail::FindCylinderAABBs())
+    .Invoke(this->CylIds,
+            this->Radii,
+            AABB.xmins,
+            AABB.ymins,
+            AABB.zmins,
+            AABB.xmaxs,
+            AABB.ymaxs,
+            AABB.zmaxs,
+            CoordsHandle);
+
+  this->SetAABBs(AABB);
+}
 
 void CylinderIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
 void CylinderIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
-template <typename Device, typename Precision>
-void CylinderIntersector::IntersectRaysImp(Device,
-                                           Ray<Precision>& rays,
-                                           bool vtkmNotUsed(returnCellIndex))
+template <typename Precision>
+void CylinderIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
 {
 
-  detail::CylinderLeafIntersector<Device> leafIntersector(this->CylIds, Radii);
+  detail::CylinderLeafWrapper leafIntersector(this->CylIds, Radii);
 
-  BVHTraverser<detail::CylinderLeafIntersector> traverser;
-  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
+  BVHTraverser traverser;
+  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
 
   RayOperations::UpdateRayStatus(rays);
 }
diff --git a/vtkm/rendering/raytracing/CylinderIntersector.h b/vtkm/rendering/raytracing/CylinderIntersector.h
index 75b1b30e2..cb7d0dddb 100644
--- a/vtkm/rendering/raytracing/CylinderIntersector.h
+++ b/vtkm/rendering/raytracing/CylinderIntersector.h
@@ -31,125 +31,6 @@ namespace raytracing
 {
 namespace detail
 {
-
-
-class FindCylinderAABBs : public vtkm::worklet::WorkletMapField
-{
-public:
-  VTKM_CONT
-  FindCylinderAABBs() {}
-  typedef void ControlSignature(FieldIn<>,
-                                FieldIn<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                WholeArrayIn<Vec3RenderingTypes>);
-  typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8, _9);
-  template <typename PointPortalType>
-  VTKM_EXEC void operator()(const vtkm::Id3 cylId,
-                            const vtkm::Float32& radius,
-                            vtkm::Float32& xmin,
-                            vtkm::Float32& ymin,
-                            vtkm::Float32& zmin,
-                            vtkm::Float32& xmax,
-                            vtkm::Float32& ymax,
-                            vtkm::Float32& zmax,
-                            const PointPortalType& points) const
-  {
-    // cast to Float32
-    vtkm::Vec<vtkm::Float32, 3> point1, point2;
-    vtkm::Vec<vtkm::Float32, 3> temp;
-
-    point1 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[1]));
-    point2 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[2]));
-
-    temp[0] = radius;
-    temp[1] = 0.0f;
-    temp[2] = 0.0f;
-    xmin = ymin = zmin = vtkm::Infinity32();
-    xmax = ymax = zmax = vtkm::NegativeInfinity32();
-
-
-    //set first point to max and min
-    Bounds(point1, radius, xmin, ymin, zmin, xmax, ymax, zmax);
-
-    Bounds(point2, radius, xmin, ymin, zmin, xmax, ymax, zmax);
-  }
-
-  VTKM_EXEC void Bounds(const vtkm::Vec<vtkm::Float32, 3>& point,
-                        const vtkm::Float32& radius,
-                        vtkm::Float32& xmin,
-                        vtkm::Float32& ymin,
-                        vtkm::Float32& zmin,
-                        vtkm::Float32& xmax,
-                        vtkm::Float32& ymax,
-                        vtkm::Float32& zmax) const
-  {
-    vtkm::Vec<vtkm::Float32, 3> temp, p;
-    temp[0] = radius;
-    temp[1] = 0.0f;
-    temp[2] = 0.0f;
-    p = point + temp;
-
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-
-    p = point - temp;
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-
-    temp[0] = 0.f;
-    temp[1] = radius;
-    temp[2] = 0.f;
-
-    p = point + temp;
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-
-    p = point - temp;
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-
-    temp[0] = 0.f;
-    temp[1] = 0.f;
-    temp[2] = radius;
-
-    p = point + temp;
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-
-    p = point - temp;
-    xmin = vtkm::Min(xmin, p[0]);
-    xmax = vtkm::Max(xmax, p[0]);
-    ymin = vtkm::Min(ymin, p[1]);
-    ymax = vtkm::Max(ymax, p[1]);
-    zmin = vtkm::Min(zmin, p[2]);
-    zmax = vtkm::Max(zmax, p[2]);
-  }
-}; //class FindAABBs
 }
 class CylinderIntersector : public ShapeIntersector
 {
@@ -164,33 +45,15 @@ public:
 
   void SetData(const vtkm::cont::CoordinateSystem& coords,
                vtkm::cont::ArrayHandle<vtkm::Id3> cylIds,
-               vtkm::cont::ArrayHandle<vtkm::Float32> radii)
-  {
-    this->Radii = radii;
-    this->CylIds = cylIds;
-    this->CoordsHandle = coords;
-    AABBs AABB;
+               vtkm::cont::ArrayHandle<vtkm::Float32> radii);
 
-    vtkm::worklet::DispatcherMapField<detail::FindCylinderAABBs>(detail::FindCylinderAABBs())
-      .Invoke(this->CylIds,
-              this->Radii,
-              AABB.xmins,
-              AABB.ymins,
-              AABB.zmins,
-              AABB.xmaxs,
-              AABB.ymaxs,
-              AABB.zmaxs,
-              CoordsHandle);
-
-    this->SetAABBs(AABB);
-  }
   void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
 
 
   void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
 
-  template <typename Device, typename Precision>
-  void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
+  template <typename Precision>
+  void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
 
 
   template <typename Precision>
diff --git a/vtkm/rendering/raytracing/MeshConnectivityBase.h b/vtkm/rendering/raytracing/MeshConnectivityBase.h
index 90685c0e0..6b41baa90 100644
--- a/vtkm/rendering/raytracing/MeshConnectivityBase.h
+++ b/vtkm/rendering/raytracing/MeshConnectivityBase.h
@@ -37,7 +37,10 @@ namespace rendering
 {
 namespace raytracing
 {
-
+//
+// Base class for different types of face-to-connecting-cell
+// and other mesh information
+//
 class VTKM_ALWAYS_EXPORT MeshConnectivityBase : public VirtualObjectBase
 {
 public:
@@ -51,6 +54,35 @@ public:
   virtual vtkm::UInt8 GetCellShape(const vtkm::Id& cellId) const = 0;
 };
 
+// A simple concrete type to wrap MeshConnectivityBase so we can
+// pass an ExeObject to worklets.
+class MeshWrapper
+{
+private:
+  MeshConnectivityBase* MeshConn;
+
+public:
+  MeshWrapper() {}
+
+  MeshWrapper(MeshConnectivityBase* meshConn)
+    : MeshConn(meshConn){};
+
+  VTKM_EXEC_CONT
+  vtkm::Id GetConnectingCell(const vtkm::Id& cellId, const vtkm::Id& face) const
+  {
+    return MeshConn->GetConnectingCell(cellId, face);
+  }
+
+  VTKM_EXEC_CONT
+  vtkm::Int32 GetCellIndices(vtkm::Id cellIndices[8], const vtkm::Id& cellId) const
+  {
+    return MeshConn->GetCellIndices(cellIndices, cellId);
+  }
+
+  VTKM_EXEC_CONT
+  vtkm::UInt8 GetCellShape(const vtkm::Id& cellId) const { return MeshConn->GetCellShape(cellId); }
+};
+
 class VTKM_ALWAYS_EXPORT MeshConnStructured : public MeshConnectivityBase
 {
 protected:
diff --git a/vtkm/rendering/raytracing/MeshConnectivityBuilder.cxx b/vtkm/rendering/raytracing/MeshConnectivityBuilder.cxx
index 9443c809e..6463c68ea 100644
--- a/vtkm/rendering/raytracing/MeshConnectivityBuilder.cxx
+++ b/vtkm/rendering/raytracing/MeshConnectivityBuilder.cxx
@@ -40,10 +40,10 @@ namespace rendering
 namespace raytracing
 {
 
-struct IsExternal
+struct IsUnique
 {
   VTKM_EXEC_CONT
-  inline bool operator()(const vtkm::Id& x) const { return (x < 0); }
+  inline bool operator()(const vtkm::Int32& x) const { return (x < 0); }
 }; //struct IsExternal
 
 class CountFaces : public vtkm::worklet::WorkletMapField
@@ -93,8 +93,9 @@ public:
                                 WholeArrayIn<>,
                                 WholeArrayIn<>,
                                 WholeArrayIn<>,
-                                WholeArrayOut<>);
-  using ExecutionSignature = void(_1, _2, WorkIndex, _3, _4, _5, _6);
+                                WholeArrayOut<>,
+                                WholeArrayInOut<>);
+  using ExecutionSignature = void(_1, _2, WorkIndex, _3, _4, _5, _6, _7);
 
   VTKM_EXEC
   inline vtkm::Int32 GetShapeOffset(const vtkm::UInt8& shapeType) const
@@ -145,14 +146,16 @@ public:
             typename ConnPortalType,
             typename ShapePortalType,
             typename OffsetPortalType,
-            typename ExternalFaceFlagType>
+            typename ExternalFaceFlagType,
+            typename UniqueFacesType>
   VTKM_EXEC inline void operator()(const MortonPortalType& mortonCodes,
                                    FaceIdPairsPortalType& faceIdPairs,
                                    const vtkm::Id& index,
                                    const ConnPortalType& connectivity,
                                    const ShapePortalType& shapes,
                                    const OffsetPortalType& offsets,
-                                   ExternalFaceFlagType& flags) const
+                                   ExternalFaceFlagType& flags,
+                                   UniqueFacesType& uniqueFaces) const
   {
     if (index == 0)
     {
@@ -252,6 +255,10 @@ public:
       flags.Set(currentIndex, myCell);
       BOUNDS_CHECK(flags, index);
       flags.Set(index, connectedCell);
+
+      // for unstructured, we want all unique faces to intersect with
+      // just choose one and mark as unique so the other gets culled.
+      uniqueFaces.Set(index, 1);
     }
   }
 }; //class Neighbor
@@ -522,7 +529,8 @@ VTKM_CONT void GenerateFaceConnnectivity(
   vtkm::cont::ArrayHandle<vtkm::Id>& faceConnectivity,
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>>& cellFaceId,
   vtkm::Float32 BoundingBox[6],
-  vtkm::cont::ArrayHandle<vtkm::Id>& faceOffsets)
+  vtkm::cont::ArrayHandle<vtkm::Id>& faceOffsets,
+  vtkm::cont::ArrayHandle<vtkm::Int32>& uniqueFaces)
 {
 
   vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
@@ -578,8 +586,8 @@ VTKM_CONT void GenerateFaceConnnectivity(
   vtkm::cont::ArrayHandle<vtkm::UInt32> faceMortonCodes;
   cellFaceId.Allocate(totalFaces);
   faceMortonCodes.Allocate(totalFaces);
-  //cellFaceId.PrepareForOutput(totalFaces, DeviceAdapter());
-  //faceMortonCodes.PrepareForOutput(totalFaces, DeviceAdapter());
+  uniqueFaces.Allocate(totalFaces);
+
   vtkm::worklet::DispatcherMapTopology<MortonCodeFace>(MortonCodeFace(inverseExtent, minPoint))
     .Invoke(cellSet, coordinates, cellOffsets, faceMortonCodes, cellFaceId);
 
@@ -594,19 +602,21 @@ VTKM_CONT void GenerateFaceConnnectivity(
   vtkm::cont::ArrayHandleConstant<vtkm::Id> negOne(-1, totalFaces);
   vtkm::cont::Algorithm::Copy(negOne, faceConnectivity);
 
+  vtkm::cont::ArrayHandleConstant<vtkm::Int32> negOne32(-1, totalFaces);
+  vtkm::cont::Algorithm::Copy(negOne32, uniqueFaces);
+
   vtkm::worklet::DispatcherMapField<MortonNeighbor>(MortonNeighbor())
-    .Invoke(faceMortonCodes, cellFaceId, conn, shapes, shapeOffsets, faceConnectivity);
+    .Invoke(faceMortonCodes, cellFaceId, conn, shapes, shapeOffsets, faceConnectivity, uniqueFaces);
 
   vtkm::Float64 time = timer.GetElapsedTime();
   Logger::GetInstance()->AddLogData("gen_face_conn", time);
 }
 
 template <typename ShapeHandleType, typename OffsetsHandleType, typename ConnHandleType>
-VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractExternalFaces(
+VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractFaces(
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>>
-    cellFaceId, // Map of cell, face, and connecting cell
-  vtkm::cont::ArrayHandle<vtkm::Id>
-    faceConnectivity, // -1 if the face does not connect to any other face
+    cellFaceId,                                     // Map of cell, face, and connecting cell
+  vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces, // -1 if the face is unique
   const ShapeHandleType& shapes,
   const ConnHandleType& conn,
   const OffsetsHandleType& shapeOffsets)
@@ -614,7 +624,7 @@ VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractExternalFaces(
 
   vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> externalFacePairs;
-  vtkm::cont::Algorithm::CopyIf(cellFaceId, faceConnectivity, externalFacePairs, IsExternal());
+  vtkm::cont::Algorithm::CopyIf(cellFaceId, uniqueFaces, externalFacePairs, IsUnique());
 
   // We need to count the number of triangle per external face
   // If it is a single cell type and it is a tet or hex, this is a special case
@@ -690,6 +700,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
 
   vtkm::cont::ArrayHandle<vtkm::Id> faceConnectivity;
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> cellFaceId;
+  vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces;
 
   GenerateFaceConnnectivity(cellSetUnstructured,
                             shapes,
@@ -699,13 +710,13 @@ void MeshConnectivityBuilder::BuildConnectivity(
                             faceConnectivity,
                             cellFaceId,
                             BoundingBox,
-                            FaceOffsets);
+                            FaceOffsets,
+                            uniqueFaces);
 
-  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> externalTriangles;
-  //External Faces
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles;
+  //Faces
 
-  externalTriangles =
-    ExtractExternalFaces(cellFaceId, faceConnectivity, shapes, conn, shapeOffsets);
+  triangles = ExtractFaces(cellFaceId, uniqueFaces, shapes, conn, shapeOffsets);
 
 
   // scatter the coonectivity into the original order
@@ -713,7 +724,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
     .Invoke(cellFaceId, this->FaceOffsets, faceConnectivity);
 
   FaceConnectivity = faceConnectivity;
-  OutsideTriangles = externalTriangles;
+  Triangles = triangles;
 
   vtkm::Float64 time = timer.GetElapsedTime();
   logger->CloseLogEntry(time);
@@ -748,6 +759,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
 
   vtkm::cont::ArrayHandle<vtkm::Id> faceConnectivity;
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> cellFaceId;
+  vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces;
 
   GenerateFaceConnnectivity(cellSetUnstructured,
                             shapes,
@@ -757,20 +769,20 @@ void MeshConnectivityBuilder::BuildConnectivity(
                             faceConnectivity,
                             cellFaceId,
                             BoundingBox,
-                            FaceOffsets);
+                            FaceOffsets,
+                            uniqueFaces);
 
-  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> externalTriangles;
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles;
   //
-  //External Faces
-  externalTriangles =
-    ExtractExternalFaces(cellFaceId, faceConnectivity, shapes, conn, shapeOffsets);
+  //Faces
+  triangles = ExtractFaces(cellFaceId, uniqueFaces, shapes, conn, shapeOffsets);
 
   // scatter the coonectivity into the original order
   vtkm::worklet::DispatcherMapField<WriteFaceConn>(WriteFaceConn())
     .Invoke(cellFaceId, this->FaceOffsets, faceConnectivity);
 
   FaceConnectivity = faceConnectivity;
-  OutsideTriangles = externalTriangles;
+  Triangles = triangles;
 
   vtkm::Float64 time = timer.GetElapsedTime();
   logger->CloseLogEntry(time);
@@ -829,9 +841,9 @@ vtkm::cont::ArrayHandle<vtkm::Id> MeshConnectivityBuilder::GetFaceOffsets()
   return FaceOffsets;
 }
 
-vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> MeshConnectivityBuilder::GetExternalTriangles()
+vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> MeshConnectivityBuilder::GetTriangles()
 {
-  return OutsideTriangles;
+  return Triangles;
 }
 
 VTKM_CONT
@@ -894,21 +906,20 @@ MeshConnContainer* MeshConnectivityBuilder::BuildConnectivity(
   {
     vtkm::cont::CellSetExplicit<> cells = cellset.Cast<vtkm::cont::CellSetExplicit<>>();
     this->BuildConnectivity(cells, coordinates.GetData(), coordBounds);
-    meshConn = new UnstructuredContainer(
-      cells, coordinates, FaceConnectivity, FaceOffsets, OutsideTriangles);
+    meshConn =
+      new UnstructuredContainer(cells, coordinates, FaceConnectivity, FaceOffsets, Triangles);
   }
   else if (type == UnstructuredSingle)
   {
     vtkm::cont::CellSetSingleType<> cells = cellset.Cast<vtkm::cont::CellSetSingleType<>>();
     this->BuildConnectivity(cells, coordinates.GetData(), coordBounds);
-    meshConn =
-      new UnstructuredSingleContainer(cells, coordinates, FaceConnectivity, OutsideTriangles);
+    meshConn = new UnstructuredSingleContainer(cells, coordinates, FaceConnectivity, Triangles);
   }
   else if (type == Structured)
   {
     vtkm::cont::CellSetStructured<3> cells = cellset.Cast<vtkm::cont::CellSetStructured<3>>();
-    OutsideTriangles = this->ExternalTrianglesStructured(cells);
-    meshConn = new StructuredContainer(cells, coordinates, OutsideTriangles);
+    Triangles = this->ExternalTrianglesStructured(cells);
+    meshConn = new StructuredContainer(cells, coordinates, Triangles);
   }
   else
   {
diff --git a/vtkm/rendering/raytracing/MeshConnectivityBuilder.h b/vtkm/rendering/raytracing/MeshConnectivityBuilder.h
index 35accd9de..7f99cc71f 100644
--- a/vtkm/rendering/raytracing/MeshConnectivityBuilder.h
+++ b/vtkm/rendering/raytracing/MeshConnectivityBuilder.h
@@ -48,7 +48,7 @@ public:
 
   vtkm::cont::ArrayHandle<vtkm::Id> GetFaceOffsets();
 
-  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetExternalTriangles();
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles();
 
 protected:
   VTKM_CONT
@@ -63,7 +63,7 @@ protected:
 
   vtkm::cont::ArrayHandle<vtkm::Id> FaceConnectivity;
   vtkm::cont::ArrayHandle<vtkm::Id> FaceOffsets;
-  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> OutsideTriangles;
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Triangles;
 };
 }
 }
diff --git a/vtkm/rendering/raytracing/MeshConnectivityContainers.cxx b/vtkm/rendering/raytracing/MeshConnectivityContainers.cxx
index 39cc29abe..c19701248 100644
--- a/vtkm/rendering/raytracing/MeshConnectivityContainers.cxx
+++ b/vtkm/rendering/raytracing/MeshConnectivityContainers.cxx
@@ -39,41 +39,28 @@ MeshConnContainer::MeshConnContainer(){};
 MeshConnContainer::~MeshConnContainer(){};
 
 template <typename T>
-VTKM_CONT void MeshConnContainer::FindEntryImpl(Ray<T>& rays,
-                                                const vtkm::cont::DeviceAdapterId deviceId)
+VTKM_CONT void MeshConnContainer::FindEntryImpl(Ray<T>& rays)
 {
   bool getCellIndex = true;
 
   Intersector.SetUseWaterTight(true);
 
-  switch (deviceId.GetValue())
-  {
-#ifdef VTKM_ENABLE_TBB
-    case VTKM_DEVICE_ADAPTER_TBB:
-      Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagTBB(), getCellIndex);
-      break;
-#endif
-#ifdef VTKM_ENABLE_CUDA
-    case VTKM_DEVICE_ADAPTER_CUDA:
-      Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagCuda(), getCellIndex);
-      break;
-#endif
-    default:
-      Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagSerial(), getCellIndex);
-      break;
-  }
+  Intersector.IntersectRays(rays, getCellIndex);
 }
 
-void MeshConnContainer::FindEntry(Ray<vtkm::Float32>& rays,
-                                  const vtkm::cont::DeviceAdapterId deviceId)
+MeshWrapper MeshConnContainer::PrepareForExecution(const vtkm::cont::DeviceAdapterId deviceId)
 {
-  this->FindEntryImpl(rays, deviceId);
+  return MeshWrapper(const_cast<MeshConnectivityBase*>(this->Construct(deviceId)));
 }
 
-void MeshConnContainer::FindEntry(Ray<vtkm::Float64>& rays,
-                                  const vtkm::cont::DeviceAdapterId deviceId)
+void MeshConnContainer::FindEntry(Ray<vtkm::Float32>& rays)
 {
-  this->FindEntryImpl(rays, deviceId);
+  this->FindEntryImpl(rays);
+}
+
+void MeshConnContainer::FindEntry(Ray<vtkm::Float64>& rays)
+{
+  this->FindEntryImpl(rays);
 }
 
 VTKM_CONT
@@ -81,13 +68,13 @@ UnstructuredContainer::UnstructuredContainer(const vtkm::cont::CellSetExplicit<>
                                              const vtkm::cont::CoordinateSystem& coords,
                                              IdHandle& faceConn,
                                              IdHandle& faceOffsets,
-                                             Id4Handle& externalTriangles)
+                                             Id4Handle& triangles)
   : FaceConnectivity(faceConn)
   , FaceOffsets(faceOffsets)
   , Cellset(cellset)
   , Coords(coords)
 {
-  this->ExternalTriangles = externalTriangles;
+  this->Triangles = triangles;
   //
   // Grab the cell arrays
   //
@@ -97,7 +84,7 @@ UnstructuredContainer::UnstructuredContainer(const vtkm::cont::CellSetExplicit<>
     Cellset.GetIndexOffsetArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
   Shapes = Cellset.GetShapesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
 
-  Intersector.SetData(Coords, ExternalTriangles);
+  Intersector.SetData(Coords, Triangles);
 }
 
 UnstructuredContainer::~UnstructuredContainer(){};
@@ -107,6 +94,19 @@ const MeshConnectivityBase* UnstructuredContainer::Construct(
 {
   switch (deviceId.GetValue())
   {
+#ifdef VTKM_ENABLE_OPENMP
+    case VTKM_DEVICE_ADAPTER_OPENMP:
+      using OMP = vtkm::cont::DeviceAdapterTagOpenMP;
+      {
+        MeshConnUnstructured<OMP> conn(this->FaceConnectivity,
+                                       this->FaceOffsets,
+                                       this->CellConn,
+                                       this->CellOffsets,
+                                       this->Shapes);
+        Handle = make_MeshConnHandle(conn);
+      }
+      return Handle.PrepareForExecution(OMP());
+#endif
 #ifdef VTKM_ENABLE_TBB
     case VTKM_DEVICE_ADAPTER_TBB:
       using TBB = vtkm::cont::DeviceAdapterTagTBB;
@@ -157,13 +157,13 @@ UnstructuredSingleContainer::UnstructuredSingleContainer(
   const vtkm::cont::CellSetSingleType<>& cellset,
   const vtkm::cont::CoordinateSystem& coords,
   IdHandle& faceConn,
-  Id4Handle& externalTriangles)
+  Id4Handle& triangles)
   : FaceConnectivity(faceConn)
   , Coords(coords)
   , Cellset(cellset)
 {
 
-  this->ExternalTriangles = externalTriangles;
+  this->Triangles = triangles;
 
   this->Intersector.SetUseWaterTight(true);
 
@@ -192,7 +192,7 @@ UnstructuredSingleContainer::UnstructuredSingleContainer(
   logger->OpenLogEntry("mesh_conn_construction");
   vtkm::cont::Timer<cont::DeviceAdapterTagSerial> timer;
 
-  Intersector.SetData(Coords, ExternalTriangles);
+  Intersector.SetData(Coords, Triangles);
 }
 
 const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
@@ -200,6 +200,20 @@ const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
 {
   switch (deviceId.GetValue())
   {
+#ifdef VTKM_ENABLE_OPENMP
+    case VTKM_DEVICE_ADAPTER_OPENMP:
+      using OMP = vtkm::cont::DeviceAdapterTagOpenMP;
+      {
+        MeshConnSingleType<OMP> conn(this->FaceConnectivity,
+                                     this->CellConnectivity,
+                                     this->CellOffsets,
+                                     this->ShapeId,
+                                     this->NumIndices,
+                                     this->NumFaces);
+        Handle = make_MeshConnHandle(conn);
+      }
+      return Handle.PrepareForExecution(OMP());
+#endif
 #ifdef VTKM_ENABLE_TBB
     case VTKM_DEVICE_ADAPTER_TBB:
       using TBB = vtkm::cont::DeviceAdapterTagTBB;
@@ -245,18 +259,18 @@ const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
 
 StructuredContainer::StructuredContainer(const vtkm::cont::CellSetStructured<3>& cellset,
                                          const vtkm::cont::CoordinateSystem& coords,
-                                         Id4Handle& externalTriangles)
+                                         Id4Handle& triangles)
   : Coords(coords)
   , Cellset(cellset)
 {
 
-  ExternalTriangles = externalTriangles;
+  Triangles = triangles;
   Intersector.SetUseWaterTight(true);
 
   PointDims = Cellset.GetPointDimensions();
   CellDims = Cellset.GetCellDimensions();
 
-  this->Intersector.SetData(Coords, ExternalTriangles);
+  this->Intersector.SetData(Coords, Triangles);
 }
 
 const MeshConnectivityBase* StructuredContainer::Construct(
@@ -268,6 +282,10 @@ const MeshConnectivityBase* StructuredContainer::Construct(
 
   switch (deviceId.GetValue())
   {
+#ifdef VTKM_ENABLE_OPENMP
+    case VTKM_DEVICE_ADAPTER_OPENMP:
+      return Handle.PrepareForExecution(vtkm::cont::DeviceAdapterTagOpenMP());
+#endif
 #ifdef VTKM_ENABLE_TBB
     case VTKM_DEVICE_ADAPTER_TBB:
       return Handle.PrepareForExecution(vtkm::cont::DeviceAdapterTagTBB());
diff --git a/vtkm/rendering/raytracing/MeshConnectivityContainers.h b/vtkm/rendering/raytracing/MeshConnectivityContainers.h
index 42e9d30c0..d33201979 100644
--- a/vtkm/rendering/raytracing/MeshConnectivityContainers.h
+++ b/vtkm/rendering/raytracing/MeshConnectivityContainers.h
@@ -31,7 +31,7 @@ namespace rendering
 namespace raytracing
 {
 
-class MeshConnContainer
+class MeshConnContainer : vtkm::cont::ExecutionObjectBase
 {
 public:
   MeshConnContainer();
@@ -39,17 +39,19 @@ public:
 
   virtual const MeshConnectivityBase* Construct(const vtkm::cont::DeviceAdapterId deviceId) = 0;
 
+  MeshWrapper PrepareForExecution(const vtkm::cont::DeviceAdapterId deviceId);
+
   template <typename T>
-  VTKM_CONT void FindEntryImpl(Ray<T>& rays, const vtkm::cont::DeviceAdapterId deviceId);
+  VTKM_CONT void FindEntryImpl(Ray<T>& rays);
 
-  void FindEntry(Ray<vtkm::Float32>& rays, const vtkm::cont::DeviceAdapterId deviceId);
+  void FindEntry(Ray<vtkm::Float32>& rays);
 
-  void FindEntry(Ray<vtkm::Float64>& rays, const vtkm::cont::DeviceAdapterId deviceId);
+  void FindEntry(Ray<vtkm::Float64>& rays);
 
 protected:
   using Id4Handle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
   // Mesh Boundary
-  Id4Handle ExternalTriangles;
+  Id4Handle Triangles;
   TriangleIntersector Intersector;
   MeshConnHandle Handle;
 };
@@ -83,7 +85,7 @@ public:
                         const vtkm::cont::CoordinateSystem& coords,
                         IdHandle& faceConn,
                         IdHandle& faceOffsets,
-                        Id4Handle& externalTriangles);
+                        Id4Handle& triangles);
 
   virtual ~UnstructuredContainer();
 
@@ -108,7 +110,7 @@ public:
   VTKM_CONT
   StructuredContainer(const vtkm::cont::CellSetStructured<3>& cellset,
                       const vtkm::cont::CoordinateSystem& coords,
-                      Id4Handle& externalTriangles);
+                      Id4Handle& triangles);
 
   const MeshConnectivityBase* Construct(const vtkm::cont::DeviceAdapterId deviceId) override;
 
diff --git a/vtkm/rendering/raytracing/QuadIntersector.cxx b/vtkm/rendering/raytracing/QuadIntersector.cxx
index f802b1a85..aac900cdb 100644
--- a/vtkm/rendering/raytracing/QuadIntersector.cxx
+++ b/vtkm/rendering/raytracing/QuadIntersector.cxx
@@ -34,17 +34,91 @@ namespace raytracing
 namespace detail
 {
 
+#define QUAD_AABB_EPSILON 1.0e-4f
+class FindQuadAABBs : public vtkm::worklet::WorkletMapField
+{
+public:
+  VTKM_CONT
+  FindQuadAABBs() {}
+  typedef void ControlSignature(FieldIn<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                FieldOut<>,
+                                WholeArrayIn<Vec3RenderingTypes>);
+  typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
+  template <typename PointPortalType>
+  VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 5> quadId,
+                            vtkm::Float32& xmin,
+                            vtkm::Float32& ymin,
+                            vtkm::Float32& zmin,
+                            vtkm::Float32& xmax,
+                            vtkm::Float32& ymax,
+                            vtkm::Float32& zmax,
+                            const PointPortalType& points) const
+  {
+    // cast to Float32
+    vtkm::Vec<vtkm::Float32, 3> q, r, s, t;
+
+    q = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[1]));
+    r = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[2]));
+    s = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[3]));
+    t = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[4]));
+
+    xmin = q[0];
+    ymin = q[1];
+    zmin = q[2];
+    xmax = xmin;
+    ymax = ymin;
+    zmax = zmin;
+    xmin = vtkm::Min(xmin, r[0]);
+    ymin = vtkm::Min(ymin, r[1]);
+    zmin = vtkm::Min(zmin, r[2]);
+    xmax = vtkm::Max(xmax, r[0]);
+    ymax = vtkm::Max(ymax, r[1]);
+    zmax = vtkm::Max(zmax, r[2]);
+    xmin = vtkm::Min(xmin, s[0]);
+    ymin = vtkm::Min(ymin, s[1]);
+    zmin = vtkm::Min(zmin, s[2]);
+    xmax = vtkm::Max(xmax, s[0]);
+    ymax = vtkm::Max(ymax, s[1]);
+    zmax = vtkm::Max(zmax, s[2]);
+    xmin = vtkm::Min(xmin, t[0]);
+    ymin = vtkm::Min(ymin, t[1]);
+    zmin = vtkm::Min(zmin, t[2]);
+    xmax = vtkm::Max(xmax, t[0]);
+    ymax = vtkm::Max(ymax, t[1]);
+    zmax = vtkm::Max(zmax, t[2]);
+
+    vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
+    const vtkm::Float32 minEpsilon = 1e-6f;
+    xEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (xmax - xmin));
+    yEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (ymax - ymin));
+    zEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (zmax - zmin));
+
+    xmin -= xEpsilon;
+    ymin -= yEpsilon;
+    zmin -= zEpsilon;
+    xmax += xEpsilon;
+    ymax += yEpsilon;
+    zmax += zEpsilon;
+  }
+
+}; //class FindAABBs
+
 template <typename Device>
-class QuadLeafIntersector : public vtkm::cont::ExecutionObjectBase
+class QuadLeafIntersector
 {
 public:
   using IdType = vtkm::Vec<vtkm::Id, 5>;
   using IdHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>>;
   using IdArrayPortal = typename IdHandle::ExecutionTypes<Device>::PortalConst;
-  using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
-  using FloatPortal = typename FloatHandle::ExecutionTypes<Device>::PortalConst;
   IdArrayPortal QuadIds;
 
+  QuadLeafIntersector() {}
+
   QuadLeafIntersector(const IdHandle& quadIds)
     : QuadIds(quadIds.PrepareForInput(Device()))
   {
@@ -228,17 +302,23 @@ public:
   }
 };
 
-struct IntersectFunctor
+class QuadExecWrapper : public vtkm::cont::ExecutionObjectBase
 {
-  template <typename Device, typename Precision>
-  VTKM_CONT bool operator()(Device,
-                            QuadIntersector* self,
-                            Ray<Precision>& rays,
-                            bool returnCellIndex)
+protected:
+  using IdType = vtkm::Vec<vtkm::Id, 5>;
+  using IdHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>>;
+  IdHandle QuadIds;
+
+public:
+  QuadExecWrapper(IdHandle& quadIds)
+    : QuadIds(quadIds)
   {
-    VTKM_IS_DEVICE_ADAPTER_TAG(Device);
-    self->IntersectRaysImp(Device(), rays, returnCellIndex);
-    return true;
+  }
+
+  template <typename Device>
+  VTKM_CONT QuadLeafIntersector<Device> PrepareForExecution(Device) const
+  {
+    return QuadLeafIntersector<Device>(QuadIds);
   }
 };
 
@@ -343,24 +423,22 @@ QuadIntersector::~QuadIntersector()
 
 void QuadIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
 void QuadIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
-template <typename Device, typename Precision>
-void QuadIntersector::IntersectRaysImp(Device,
-                                       Ray<Precision>& rays,
-                                       bool vtkmNotUsed(returnCellIndex))
+template <typename Precision>
+void QuadIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
 {
 
-  detail::QuadLeafIntersector<Device> leafIntersector(this->QuadIds);
+  detail::QuadExecWrapper leafIntersector(this->QuadIds);
 
-  BVHTraverser<detail::QuadLeafIntersector> traverser;
-  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
+  BVHTraverser traverser;
+  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
 
   RayOperations::UpdateRayStatus(rays);
 }
@@ -401,6 +479,27 @@ void QuadIntersector::IntersectionData(Ray<vtkm::Float64>& rays,
   IntersectionDataImp(rays, scalarField, scalarRange);
 }
 
+void QuadIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
+                              vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds)
+{
+
+  this->QuadIds = quadIds;
+  this->CoordsHandle = coords;
+  AABBs AABB;
+
+  vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker;
+  faabbsInvoker.Invoke(this->QuadIds,
+                       AABB.xmins,
+                       AABB.ymins,
+                       AABB.zmins,
+                       AABB.xmaxs,
+                       AABB.ymaxs,
+                       AABB.zmaxs,
+                       CoordsHandle);
+
+  this->SetAABBs(AABB);
+}
+
 vtkm::Id QuadIntersector::GetNumberOfShapes() const
 {
   return QuadIds.GetNumberOfValues();
diff --git a/vtkm/rendering/raytracing/QuadIntersector.h b/vtkm/rendering/raytracing/QuadIntersector.h
index fbfb1dd2d..091baef4e 100644
--- a/vtkm/rendering/raytracing/QuadIntersector.h
+++ b/vtkm/rendering/raytracing/QuadIntersector.h
@@ -31,80 +31,6 @@ namespace raytracing
 {
 namespace detail
 {
-
-#define AABB_EPSILON 1.0e-4f
-class FindQuadAABBs : public vtkm::worklet::WorkletMapField
-{
-public:
-  VTKM_CONT
-  FindQuadAABBs() {}
-  typedef void ControlSignature(FieldIn<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                WholeArrayIn<Vec3RenderingTypes>);
-  typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
-  template <typename PointPortalType>
-  VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 5> quadId,
-                            vtkm::Float32& xmin,
-                            vtkm::Float32& ymin,
-                            vtkm::Float32& zmin,
-                            vtkm::Float32& xmax,
-                            vtkm::Float32& ymax,
-                            vtkm::Float32& zmax,
-                            const PointPortalType& points) const
-  {
-    // cast to Float32
-    vtkm::Vec<vtkm::Float32, 3> q, r, s, t;
-
-    q = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[1]));
-    r = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[2]));
-    s = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[3]));
-    t = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[4]));
-
-    xmin = q[0];
-    ymin = q[1];
-    zmin = q[2];
-    xmax = xmin;
-    ymax = ymin;
-    zmax = zmin;
-    xmin = vtkm::Min(xmin, r[0]);
-    ymin = vtkm::Min(ymin, r[1]);
-    zmin = vtkm::Min(zmin, r[2]);
-    xmax = vtkm::Max(xmax, r[0]);
-    ymax = vtkm::Max(ymax, r[1]);
-    zmax = vtkm::Max(zmax, r[2]);
-    xmin = vtkm::Min(xmin, s[0]);
-    ymin = vtkm::Min(ymin, s[1]);
-    zmin = vtkm::Min(zmin, s[2]);
-    xmax = vtkm::Max(xmax, s[0]);
-    ymax = vtkm::Max(ymax, s[1]);
-    zmax = vtkm::Max(zmax, s[2]);
-    xmin = vtkm::Min(xmin, t[0]);
-    ymin = vtkm::Min(ymin, t[1]);
-    zmin = vtkm::Min(zmin, t[2]);
-    xmax = vtkm::Max(xmax, t[0]);
-    ymax = vtkm::Max(ymax, t[1]);
-    zmax = vtkm::Max(zmax, t[2]);
-
-    vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
-    const vtkm::Float32 minEpsilon = 1e-6f;
-    xEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (xmax - xmin));
-    yEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (ymax - ymin));
-    zEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (zmax - zmin));
-
-    xmin -= xEpsilon;
-    ymin -= yEpsilon;
-    zmin -= zEpsilon;
-    xmax += xEpsilon;
-    ymax += yEpsilon;
-    zmax += zEpsilon;
-  }
-
-}; //class FindAABBs
 }
 class QuadIntersector : public ShapeIntersector
 {
@@ -117,42 +43,15 @@ public:
 
 
   void SetData(const vtkm::cont::CoordinateSystem& coords,
-               vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds)
-  {
+               vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds);
 
-    this->QuadIds = quadIds;
-    this->CoordsHandle = coords;
-    AABBs AABB;
-
-    vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker;
-    faabbsInvoker.Invoke(this->QuadIds,
-                         AABB.xmins,
-                         AABB.ymins,
-                         AABB.zmins,
-                         AABB.xmaxs,
-                         AABB.ymaxs,
-                         AABB.zmaxs,
-                         CoordsHandle);
-
-    //vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker(detail::FindQuadAABBs())
-    //  .Invoke(this->QuadIds,
-    //          AABB.xmins,
-    //          AABB.ymins,
-    //          AABB.zmins,
-    //          AABB.xmaxs,
-    //          AABB.ymaxs,
-    //          AABB.zmaxs,
-    //          CoordsHandle);
-
-    this->SetAABBs(AABB);
-  }
   void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
 
 
   void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
 
-  template <typename Device, typename Precision>
-  void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
+  template <typename Precision>
+  void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
 
 
   template <typename Precision>
diff --git a/vtkm/rendering/raytracing/RayOperations.h b/vtkm/rendering/raytracing/RayOperations.h
index 434e5aba8..0637f1dc4 100644
--- a/vtkm/rendering/raytracing/RayOperations.h
+++ b/vtkm/rendering/raytracing/RayOperations.h
@@ -113,11 +113,11 @@ public:
 class RayOperations
 {
 public:
-  template <typename Device, typename T>
-  static void ResetStatus(Ray<T>& rays, vtkm::UInt8 status, Device)
+  template <typename T>
+  static void ResetStatus(Ray<T>& rays, vtkm::UInt8 status)
   {
     vtkm::cont::ArrayHandleConstant<vtkm::UInt8> statusHandle(status, rays.NumRays);
-    vtkm::cont::Algorithm::Copy(Device(), statusHandle, rays.Status);
+    vtkm::cont::Algorithm::Copy(statusHandle, rays.Status);
   }
 
   //
@@ -146,8 +146,8 @@ public:
                               const vtkm::rendering::Camera& camera,
                               const vtkm::rendering::CanvasRayTracer& canvas);
 
-  template <typename Device, typename T>
-  static vtkm::Id RaysInMesh(Ray<T>& rays, Device)
+  template <typename T>
+  static vtkm::Id RaysInMesh(Ray<T>& rays)
   {
     vtkm::Vec<UInt8, 2> maskValues;
     maskValues[0] = RAY_ACTIVE;
@@ -157,17 +157,16 @@ public:
 
     vtkm::worklet::DispatcherMapField<ManyMask<vtkm::UInt8, 2>> dispatcher{ (
       ManyMask<vtkm::UInt8, 2>{ maskValues }) };
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(rays.Status, masks);
     vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
     const vtkm::Id initVal = 0;
-    vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
+    vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
 
     return count;
   }
 
-  template <typename Device, typename T>
-  static vtkm::Id GetStatusCount(Ray<T>& rays, vtkm::Id status, Device)
+  template <typename T>
+  static vtkm::Id GetStatusCount(Ray<T>& rays, vtkm::Id status)
   {
     vtkm::UInt8 statusUInt8;
     if (status < 0 || status > 255)
@@ -180,17 +179,16 @@ public:
 
     vtkm::worklet::DispatcherMapField<Mask<vtkm::UInt8>> dispatcher{ (
       Mask<vtkm::UInt8>{ statusUInt8 }) };
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(rays.Status, masks);
     vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
     const vtkm::Id initVal = 0;
-    vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
+    vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
 
     return count;
   }
 
-  template <typename Device, typename T>
-  static vtkm::Id RaysProcessed(Ray<T>& rays, Device)
+  template <typename T>
+  static vtkm::Id RaysProcessed(Ray<T>& rays)
   {
     vtkm::Vec<UInt8, 3> maskValues;
     maskValues[0] = RAY_TERMINATED;
@@ -201,17 +199,16 @@ public:
 
     vtkm::worklet::DispatcherMapField<ManyMask<vtkm::UInt8, 3>> dispatcher{ (
       ManyMask<vtkm::UInt8, 3>{ maskValues }) };
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(rays.Status, masks);
     vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
     const vtkm::Id initVal = 0;
-    vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
+    vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
 
     return count;
   }
 
-  template <typename Device, typename T>
-  static vtkm::cont::ArrayHandle<vtkm::UInt8> CompactActiveRays(Ray<T>& rays, Device)
+  template <typename T>
+  static vtkm::cont::ArrayHandle<vtkm::UInt8> CompactActiveRays(Ray<T>& rays)
   {
     vtkm::Vec<UInt8, 1> maskValues;
     maskValues[0] = RAY_ACTIVE;
@@ -220,7 +217,6 @@ public:
 
     vtkm::worklet::DispatcherMapField<Mask<vtkm::UInt8>> dispatcher{ (
       Mask<vtkm::UInt8>{ statusUInt8 }) };
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(rays.Status, masks);
 
     vtkm::cont::ArrayHandle<T> emptyHandle;
@@ -261,7 +257,7 @@ public:
         break;
       }
       vtkm::cont::ArrayHandle<T> compacted;
-      vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(*floatArrayPointers[i], masks, compacted);
+      vtkm::cont::Algorithm::CopyIf(*floatArrayPointers[i], masks, compacted);
       *floatArrayPointers[i] = compacted;
     }
 
@@ -275,15 +271,15 @@ public:
     rays.Dir = vtkm::cont::make_ArrayHandleCompositeVector(rays.DirX, rays.DirY, rays.DirZ);
 
     vtkm::cont::ArrayHandle<vtkm::Id> compactedHits;
-    vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.HitIdx, masks, compactedHits);
+    vtkm::cont::Algorithm::CopyIf(rays.HitIdx, masks, compactedHits);
     rays.HitIdx = compactedHits;
 
     vtkm::cont::ArrayHandle<vtkm::Id> compactedPixels;
-    vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.PixelIdx, masks, compactedPixels);
+    vtkm::cont::Algorithm::CopyIf(rays.PixelIdx, masks, compactedPixels);
     rays.PixelIdx = compactedPixels;
 
     vtkm::cont::ArrayHandle<vtkm::UInt8> compactedStatus;
-    vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.Status, masks, compactedStatus);
+    vtkm::cont::Algorithm::CopyIf(rays.Status, masks, compactedStatus);
     rays.Status = compactedStatus;
 
     rays.NumRays = rays.Status.GetPortalConstControl().GetNumberOfValues();
@@ -291,7 +287,7 @@ public:
     const size_t bufferCount = static_cast<size_t>(rays.Buffers.size());
     for (size_t i = 0; i < bufferCount; ++i)
     {
-      ChannelBufferOperations::Compact(rays.Buffers[i], masks, rays.NumRays, Device());
+      ChannelBufferOperations::Compact(rays.Buffers[i], masks, rays.NumRays);
     }
     return masks;
   }
@@ -340,12 +336,11 @@ public:
     }
   }
 
-  template <typename Device, typename T>
-  static void CopyDistancesToMin(Ray<T> rays, Device, const T offset = 0.f)
+  template <typename T>
+  static void CopyDistancesToMin(Ray<T> rays, const T offset = 0.f)
   {
     vtkm::worklet::DispatcherMapField<CopyAndOffsetMask<T>> dispatcher{ (
       CopyAndOffsetMask<T>{ offset, RAY_EXITED_MESH }) };
-    dispatcher.SetDevice(Device());
     dispatcher.Invoke(rays.Distance, rays.MinDistance, rays.Status);
   }
 };
diff --git a/vtkm/rendering/raytracing/SphereIntersector.cxx b/vtkm/rendering/raytracing/SphereIntersector.cxx
index af587af37..22bc05d1f 100644
--- a/vtkm/rendering/raytracing/SphereIntersector.cxx
+++ b/vtkm/rendering/raytracing/SphereIntersector.cxx
@@ -141,6 +141,8 @@ public:
   IdArrayPortal PointIds;
   FloatPortal Radii;
 
+  SphereLeafIntersector() {}
+
   SphereLeafIntersector(const IdHandle& pointIds, const FloatHandle& radii)
     : PointIds(pointIds.PrepareForInput(Device()))
     , Radii(radii.PrepareForInput(Device()))
@@ -193,17 +195,25 @@ public:
   }
 };
 
-struct IntersectFunctor
+class SphereLeafWrapper : public vtkm::cont::ExecutionObjectBase
 {
-  template <typename Device, typename Precision>
-  VTKM_CONT bool operator()(Device,
-                            SphereIntersector* self,
-                            Ray<Precision>& rays,
-                            bool returnCellIndex)
+protected:
+  using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
+  using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
+  IdHandle PointIds;
+  FloatHandle Radii;
+
+public:
+  SphereLeafWrapper(IdHandle& pointIds, FloatHandle radii)
+    : PointIds(pointIds)
+    , Radii(radii)
   {
-    VTKM_IS_DEVICE_ADAPTER_TAG(Device);
-    self->IntersectRaysImp(Device(), rays, returnCellIndex);
-    return true;
+  }
+
+  template <typename Device>
+  VTKM_CONT SphereLeafIntersector<Device> PrepareForExecution(Device) const
+  {
+    return SphereLeafIntersector<Device>(PointIds, Radii);
   }
 };
 
@@ -321,24 +331,22 @@ void SphereIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
 
 void SphereIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
 void SphereIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
 {
-  vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
+  IntersectRaysImp(rays, returnCellIndex);
 }
 
-template <typename Device, typename Precision>
-void SphereIntersector::IntersectRaysImp(Device,
-                                         Ray<Precision>& rays,
-                                         bool vtkmNotUsed(returnCellIndex))
+template <typename Precision>
+void SphereIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
 {
 
-  detail::SphereLeafIntersector<Device> leafIntersector(this->PointIds, Radii);
+  detail::SphereLeafWrapper leafIntersector(this->PointIds, Radii);
 
-  BVHTraverser<detail::SphereLeafIntersector> traverser;
-  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
+  BVHTraverser traverser;
+  traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
 
   RayOperations::UpdateRayStatus(rays);
 }
diff --git a/vtkm/rendering/raytracing/SphereIntersector.h b/vtkm/rendering/raytracing/SphereIntersector.h
index d94715510..742ea90ad 100644
--- a/vtkm/rendering/raytracing/SphereIntersector.h
+++ b/vtkm/rendering/raytracing/SphereIntersector.h
@@ -48,8 +48,8 @@ public:
 
   void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
 
-  template <typename Device, typename Precision>
-  void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
+  template <typename Precision>
+  void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
 
 
   template <typename Precision>
diff --git a/vtkm/rendering/raytracing/TriangleIntersector.h b/vtkm/rendering/raytracing/TriangleIntersector.h
index 67ef2a40d..954ce2d57 100644
--- a/vtkm/rendering/raytracing/TriangleIntersector.h
+++ b/vtkm/rendering/raytracing/TriangleIntersector.h
@@ -19,16 +19,11 @@
 //============================================================================
 #ifndef vtk_m_rendering_raytracing_TriagnleIntersector_h
 #define vtk_m_rendering_raytracing_TriagnleIntersector_h
-#include <cstring>
-#include <vtkm/cont/ArrayHandle.h>
-#include <vtkm/cont/ArrayHandleCompositeVector.h>
-#include <vtkm/rendering/internal/RunTriangulator.h>
-#include <vtkm/rendering/raytracing/BVHTraverser.h>
-#include <vtkm/rendering/raytracing/BoundingVolumeHierarchy.h>
+
+#include <vtkm/cont/DataSet.h>
 #include <vtkm/rendering/raytracing/Ray.h>
-#include <vtkm/rendering/raytracing/RayOperations.h>
 #include <vtkm/rendering/raytracing/ShapeIntersector.h>
-#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/rendering/vtkm_rendering_export.h>
 
 namespace vtkm
 {
@@ -37,896 +32,44 @@ namespace rendering
 namespace raytracing
 {
 
-class Moller
-{
-public:
-  template <typename Precision>
-  VTKM_EXEC void IntersectTri(const vtkm::Vec<Precision, 3>& a,
-                              const vtkm::Vec<Precision, 3>& b,
-                              const vtkm::Vec<Precision, 3>& c,
-                              const vtkm::Vec<Precision, 3>& dir,
-                              Precision& distance,
-                              Precision& u,
-                              Precision& v,
-                              const vtkm::Vec<Precision, 3>& origin) const
-  {
-    const vtkm::Float32 EPSILON2 = 0.0001f;
-
-    vtkm::Vec<Precision, 3> e1 = b - a;
-    vtkm::Vec<Precision, 3> e2 = c - a;
-
-    vtkm::Vec<Precision, 3> p;
-    p[0] = dir[1] * e2[2] - dir[2] * e2[1];
-    p[1] = dir[2] * e2[0] - dir[0] * e2[2];
-    p[2] = dir[0] * e2[1] - dir[1] * e2[0];
-    Precision dot = e1[0] * p[0] + e1[1] * p[1] + e1[2] * p[2];
-    if (dot != 0.f)
-    {
-      dot = 1.f / dot;
-      vtkm::Vec<Precision, 3> t;
-      t = origin - a;
-
-      u = (t[0] * p[0] + t[1] * p[1] + t[2] * p[2]) * dot;
-      if (u >= (0.f - EPSILON2) && u <= (1.f + EPSILON2))
-      {
-
-        vtkm::Vec<Precision, 3> q; // = t % e1;
-        q[0] = t[1] * e1[2] - t[2] * e1[1];
-        q[1] = t[2] * e1[0] - t[0] * e1[2];
-        q[2] = t[0] * e1[1] - t[1] * e1[0];
-        v = (dir[0] * q[0] + dir[1] * q[1] + dir[2] * q[2]) * dot;
-        if (v >= (0.f - EPSILON2) && v <= (1.f + EPSILON2) && !(u + v > 1.f))
-        {
-          distance = (e2[0] * q[0] + e2[1] * q[1] + e2[2] * q[2]) * dot;
-        }
-      }
-    }
-  }
-
-}; //Moller
-
-
-
-// TODO: optimization for sorting ray dims before this call.
-//       This is called multiple times and kz,kx, and ky are
-//       constant for the ray
-
-
-class WaterTight
-{
-public:
-  template <typename Precision>
-  VTKM_EXEC inline void FindDir(const vtkm::Vec<Precision, 3>& dir,
-                                vtkm::Vec<Precision, 3>& s,
-                                vtkm::Vec<Int32, 3>& k) const
-  {
-    //Find max ray direction
-    k[2] = 0;
-    if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
-    {
-      if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
-        k[2] = 0;
-      else
-        k[2] = 2;
-    }
-    else
-    {
-      if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
-        k[2] = 1;
-      else
-        k[2] = 2;
-    }
-
-    k[0] = k[2] + 1;
-    if (k[0] == 3)
-      k[0] = 0;
-    k[1] = k[0] + 1;
-    if (k[1] == 3)
-      k[1] = 0;
-
-    if (dir[k[2]] < 0.f)
-    {
-      vtkm::Int32 temp = k[1];
-      k[1] = k[0];
-      k[0] = temp;
-    }
-
-    s[0] = dir[k[0]] / dir[k[2]];
-    s[1] = dir[k[1]] / dir[k[2]];
-    s[2] = 1.f / dir[k[2]];
-  }
-
-  template <typename Precision>
-  VTKM_EXEC_CONT inline void IntersectTri(const vtkm::Vec<Precision, 3>& a,
-                                          const vtkm::Vec<Precision, 3>& b,
-                                          const vtkm::Vec<Precision, 3>& c,
-                                          const vtkm::Vec<Precision, 3>& dir,
-                                          Precision& distance,
-                                          Precision& u,
-                                          Precision& v,
-                                          const vtkm::Vec<Precision, 3>& origin) const
-  {
-    vtkm::Vec<Int32, 3> k;
-    vtkm::Vec<Precision, 3> s;
-    //Find max ray direction
-    k[2] = 0;
-    if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
-    {
-      if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
-        k[2] = 0;
-      else
-        k[2] = 2;
-    }
-    else
-    {
-      if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
-        k[2] = 1;
-      else
-        k[2] = 2;
-    }
-
-    k[0] = k[2] + 1;
-    if (k[0] == 3)
-      k[0] = 0;
-    k[1] = k[0] + 1;
-    if (k[1] == 3)
-      k[1] = 0;
-
-    if (dir[k[2]] < 0.f)
-    {
-      vtkm::Int32 temp = k[1];
-      k[1] = k[0];
-      k[0] = temp;
-    }
-
-    s[0] = dir[k[0]] / dir[k[2]];
-    s[1] = dir[k[1]] / dir[k[2]];
-    s[2] = 1.f / dir[k[2]];
-
-    vtkm::Vec<Precision, 3> A, B, C;
-    A = a - origin;
-    B = b - origin;
-    C = c - origin;
-
-    const Precision Ax = A[k[0]] - s[0] * A[k[2]];
-    const Precision Ay = A[k[1]] - s[1] * A[k[2]];
-    const Precision Bx = B[k[0]] - s[0] * B[k[2]];
-    const Precision By = B[k[1]] - s[1] * B[k[2]];
-    const Precision Cx = C[k[0]] - s[0] * C[k[2]];
-    const Precision Cy = C[k[1]] - s[1] * C[k[2]];
-
-    //scaled barycentric coords
-    u = Cx * By - Cy * Bx;
-    v = Ax * Cy - Ay * Cx;
-    Precision w = Bx * Ay - By * Ax;
-    if (u == 0.f || v == 0.f || w == 0.f)
-    {
-      vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
-      vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
-      u = vtkm::Float32(CxBy - CyBx);
-
-      vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
-      vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
-      v = vtkm::Float32(AxCy - AyCx);
-
-      vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
-      vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
-      w = vtkm::Float32(BxAy - ByAx);
-    }
-    Precision low = vtkm::Min(u, vtkm::Min(v, w));
-    Precision high = vtkm::Max(u, vtkm::Max(v, w));
-
-    bool invalid = (low < 0.) && (high > 0.);
-
-    Precision det = u + v + w;
-
-    if (det == 0.)
-      invalid = true;
-
-    const Precision Az = s[2] * A[k[2]];
-    const Precision Bz = s[2] * B[k[2]];
-    const Precision Cz = s[2] * C[k[2]];
-
-    det = 1.f / det;
-
-    u = u * det;
-    v = v * det;
-
-    distance = (u * Az + v * Bz + w * det * Cz);
-    u = v;
-    v = w * det;
-    if (invalid)
-      distance = -1.;
-  }
-
-  template <typename Precision>
-  VTKM_EXEC inline void IntersectTriSn(const vtkm::Vec<Precision, 3>& a,
-                                       const vtkm::Vec<Precision, 3>& b,
-                                       const vtkm::Vec<Precision, 3>& c,
-                                       const vtkm::Vec<Precision, 3>& s,
-                                       const vtkm::Vec<Int32, 3>& k,
-                                       Precision& distance,
-                                       Precision& u,
-                                       Precision& v,
-                                       const vtkm::Vec<Precision, 3>& origin) const
-  {
-    vtkm::Vec<Precision, 3> A, B, C;
-    A = a - origin;
-    B = b - origin;
-    C = c - origin;
-
-    const Precision Ax = A[k[0]] - s[0] * A[k[2]];
-    const Precision Ay = A[k[1]] - s[1] * A[k[2]];
-    const Precision Bx = B[k[0]] - s[0] * B[k[2]];
-    const Precision By = B[k[1]] - s[1] * B[k[2]];
-    const Precision Cx = C[k[0]] - s[0] * C[k[2]];
-    const Precision Cy = C[k[1]] - s[1] * C[k[2]];
-
-    //scaled barycentric coords
-    u = Cx * By - Cy * Bx;
-    v = Ax * Cy - Ay * Cx;
-    Precision w = Bx * Ay - By * Ax;
-    if (u == 0.f || v == 0.f || w == 0.f)
-    {
-      vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
-      vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
-      u = vtkm::Float32(CxBy - CyBx);
-
-      vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
-      vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
-      v = vtkm::Float32(AxCy - AyCx);
-
-      vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
-      vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
-      w = vtkm::Float32(BxAy - ByAx);
-    }
-
-    Precision low = vtkm::Min(u, vtkm::Min(v, w));
-    Precision high = vtkm::Max(u, vtkm::Max(v, w));
-
-    bool invalid = (low < 0.) && (high > 0.);
-
-    Precision det = u + v + w;
-
-    if (det == 0.)
-      invalid = true;
-
-    const Precision Az = s[2] * A[k[2]];
-    const Precision Bz = s[2] * B[k[2]];
-    const Precision Cz = s[2] * C[k[2]];
-
-    det = 1.f / det;
-
-    u = u * det;
-    v = v * det;
-
-    distance = (u * Az + v * Bz + w * det * Cz);
-    u = v;
-    v = w * det;
-    if (invalid)
-      distance = -1.;
-  }
-}; //WaterTight
-
-template <>
-VTKM_EXEC inline void WaterTight::IntersectTri<vtkm::Float64>(
-  const vtkm::Vec<vtkm::Float64, 3>& a,
-  const vtkm::Vec<vtkm::Float64, 3>& b,
-  const vtkm::Vec<vtkm::Float64, 3>& c,
-  const vtkm::Vec<vtkm::Float64, 3>& dir,
-  vtkm::Float64& distance,
-  vtkm::Float64& u,
-  vtkm::Float64& v,
-  const vtkm::Vec<vtkm::Float64, 3>& origin) const
-{
-  //Find max ray direction
-  int kz = 0;
-  if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
-  {
-    if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
-      kz = 0;
-    else
-      kz = 2;
-  }
-  else
-  {
-    if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
-      kz = 1;
-    else
-      kz = 2;
-  }
-
-  vtkm::Int32 kx = kz + 1;
-  if (kx == 3)
-    kx = 0;
-  vtkm::Int32 ky = kx + 1;
-  if (ky == 3)
-    ky = 0;
-
-  if (dir[kz] < 0.f)
-  {
-    vtkm::Int32 temp = ky;
-    ky = kx;
-    kx = temp;
-  }
-
-  vtkm::Float64 Sx = dir[kx] / dir[kz];
-  vtkm::Float64 Sy = dir[ky] / dir[kz];
-  vtkm::Float64 Sz = 1. / dir[kz];
-
-
-
-  vtkm::Vec<vtkm::Float64, 3> A, B, C;
-  A = a - origin;
-  B = b - origin;
-  C = c - origin;
-
-  const vtkm::Float64 Ax = A[kx] - Sx * A[kz];
-  const vtkm::Float64 Ay = A[ky] - Sy * A[kz];
-  const vtkm::Float64 Bx = B[kx] - Sx * B[kz];
-  const vtkm::Float64 By = B[ky] - Sy * B[kz];
-  const vtkm::Float64 Cx = C[kx] - Sx * C[kz];
-  const vtkm::Float64 Cy = C[ky] - Sy * C[kz];
-
-  //scaled barycentric coords
-  u = Cx * By - Cy * Bx;
-  v = Ax * Cy - Ay * Cx;
-
-  vtkm::Float64 w = Bx * Ay - By * Ax;
-
-  vtkm::Float64 low = vtkm::Min(u, vtkm::Min(v, w));
-  vtkm::Float64 high = vtkm::Max(u, vtkm::Max(v, w));
-  bool invalid = (low < 0.) && (high > 0.);
-
-  vtkm::Float64 det = u + v + w;
-
-  if (det == 0.)
-    invalid = true;
-
-  const vtkm::Float64 Az = Sz * A[kz];
-  const vtkm::Float64 Bz = Sz * B[kz];
-  const vtkm::Float64 Cz = Sz * C[kz];
-
-  det = 1. / det;
-
-  u = u * det;
-  v = v * det;
-
-  distance = (u * Az + v * Bz + w * det * Cz);
-  u = v;
-  v = w * det;
-  if (invalid)
-    distance = -1.;
-}
-
-template <typename Device>
-//class WaterTightTriLeafIntersector : public vtkm::exec::ExecutionObjectBase
-class WaterTightTriLeafIntersector
-{
-public:
-  using Id4Handle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
-  using Id4ArrayPortal = typename Id4Handle::ExecutionTypes<Device>::PortalConst;
-  Id4ArrayPortal Triangles;
-
-public:
-  WaterTightTriLeafIntersector(const Id4Handle& triangles)
-    : Triangles(triangles.PrepareForInput(Device()))
-  {
-  }
-
-  WaterTightTriLeafIntersector(const WaterTightTriLeafIntersector<Device>& other)
-    : Triangles(other.Triangles)
-  {
-  }
-
-  template <typename PointPortalType, typename LeafPortalType, typename Precision>
-  VTKM_EXEC inline void IntersectLeaf(const vtkm::Int32& currentNode,
-                                      const vtkm::Vec<Precision, 3>& origin,
-                                      const vtkm::Vec<Precision, 3>& dir,
-                                      const PointPortalType& points,
-                                      vtkm::Id& hitIndex,
-                                      Precision& closestDistance,
-                                      Precision& minU,
-                                      Precision& minV,
-                                      LeafPortalType leafs,
-                                      const Precision& minDistance) const
-  {
-    const vtkm::Id triangleCount = leafs.Get(currentNode);
-    WaterTight intersector;
-    for (vtkm::Id i = 1; i <= triangleCount; ++i)
-    {
-      const vtkm::Id triIndex = leafs.Get(currentNode + i);
-      vtkm::Vec<Id, 4> triangle = Triangles.Get(triIndex);
-      vtkm::Vec<Precision, 3> a = vtkm::Vec<Precision, 3>(points.Get(triangle[1]));
-      vtkm::Vec<Precision, 3> b = vtkm::Vec<Precision, 3>(points.Get(triangle[2]));
-      vtkm::Vec<Precision, 3> c = vtkm::Vec<Precision, 3>(points.Get(triangle[3]));
-      Precision distance = -1.;
-      Precision u, v;
-
-      intersector.IntersectTri(a, b, c, dir, distance, u, v, origin);
-      if (distance != -1. && distance < closestDistance && distance > minDistance)
-      {
-        closestDistance = distance;
-        minU = u;
-        minV = v;
-        hitIndex = triIndex;
-      }
-    } // for
-  }
-};
-
-template <typename Device>
-//class MollerTriLeafIntersector : public vtkm::exec::ExecutionObjectBase
-class MollerTriLeafIntersector
-{
-  //protected:
-public:
-  using Id4Handle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
-  using Id4ArrayPortal = typename Id4Handle::ExecutionTypes<Device>::PortalConst;
-  Id4ArrayPortal Triangles;
-
-public:
-  MollerTriLeafIntersector(const Id4Handle& triangles)
-    : Triangles(triangles.PrepareForInput(Device()))
-  {
-  }
-
-  template <typename PointPortalType, typename LeafPortalType, typename Precision>
-  VTKM_EXEC inline void IntersectLeaf(const vtkm::Int32& currentNode,
-                                      const vtkm::Vec<Precision, 3>& origin,
-                                      const vtkm::Vec<Precision, 3>& dir,
-                                      const PointPortalType& points,
-                                      vtkm::Id& hitIndex,
-                                      Precision& closestDistance,
-                                      Precision& minU,
-                                      Precision& minV,
-                                      LeafPortalType leafs,
-                                      const Precision& minDistance) const
-  {
-    const vtkm::Id triangleCount = leafs.Get(currentNode);
-    Moller intersector;
-    for (vtkm::Id i = 1; i <= triangleCount; ++i)
-    {
-      const vtkm::Id triIndex = leafs.Get(currentNode + i);
-      vtkm::Vec<Id, 4> triangle = Triangles.Get(triIndex);
-      vtkm::Vec<Precision, 3> a = vtkm::Vec<Precision, 3>(points.Get(triangle[1]));
-      vtkm::Vec<Precision, 3> b = vtkm::Vec<Precision, 3>(points.Get(triangle[2]));
-      vtkm::Vec<Precision, 3> c = vtkm::Vec<Precision, 3>(points.Get(triangle[3]));
-      Precision distance = -1.;
-      Precision u, v;
-
-      intersector.IntersectTri(a, b, c, dir, distance, u, v, origin);
-
-      if (distance != -1. && distance < closestDistance && distance > minDistance)
-      {
-        closestDistance = distance;
-        minU = u;
-        minV = v;
-        hitIndex = triIndex;
-      }
-    } // for
-  }
-};
-
-template <typename T>
-VTKM_EXEC inline void swap(T& a, T& b)
-{
-  T tmp = a;
-  a = b;
-  b = tmp;
-}
-
-VTKM_EXEC
-inline vtkm::Float32 up(const vtkm::Float32& a)
-{
-  return (a > 0.f) ? a * (1.f + vtkm::Float32(2e-23)) : a * (1.f - vtkm::Float32(2e-23));
-}
-
-VTKM_EXEC
-inline vtkm::Float32 down(const vtkm::Float32& a)
-{
-  return (a > 0.f) ? a * (1.f - vtkm::Float32(2e-23)) : a * (1.f + vtkm::Float32(2e-23));
-}
-
-VTKM_EXEC
-inline vtkm::Float32 upFast(const vtkm::Float32& a)
-{
-  return a * (1.f + vtkm::Float32(2e-23));
-}
-
-VTKM_EXEC
-inline vtkm::Float32 downFast(const vtkm::Float32& a)
-{
-  return a * (1.f - vtkm::Float32(2e-23));
-}
-
-namespace detail
-{
-
-class TriangleIntersectionData
-{
-public:
-  // Worklet to calutate the normals of a triagle if
-  // none are stored in the data set
-  class CalculateNormals : public vtkm::worklet::WorkletMapField
-  {
-  public:
-    VTKM_CONT
-    CalculateNormals() {}
-    typedef void ControlSignature(FieldIn<>,
-                                  FieldIn<>,
-                                  FieldOut<>,
-                                  FieldOut<>,
-                                  FieldOut<>,
-                                  WholeArrayIn<Vec3RenderingTypes>,
-                                  WholeArrayIn<>);
-    typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7);
-    template <typename Precision, typename PointPortalType, typename IndicesPortalType>
-    VTKM_EXEC inline void operator()(const vtkm::Id& hitIndex,
-                                     const vtkm::Vec<Precision, 3>& rayDir,
-                                     Precision& normalX,
-                                     Precision& normalY,
-                                     Precision& normalZ,
-                                     const PointPortalType& points,
-                                     const IndicesPortalType& indicesPortal) const
-    {
-      if (hitIndex < 0)
-        return;
-
-      vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
-      vtkm::Vec<Precision, 3> a = points.Get(indices[1]);
-      vtkm::Vec<Precision, 3> b = points.Get(indices[2]);
-      vtkm::Vec<Precision, 3> c = points.Get(indices[3]);
-
-      vtkm::Vec<Precision, 3> normal = vtkm::TriangleNormal(a, b, c);
-      vtkm::Normalize(normal);
-
-      //flip the normal if its pointing the wrong way
-      if (vtkm::dot(normal, rayDir) > 0.f)
-        normal = -normal;
-      normalX = normal[0];
-      normalY = normal[1];
-      normalZ = normal[2];
-    }
-  }; //class CalculateNormals
-
-  template <typename Precision>
-  class LerpScalar : public vtkm::worklet::WorkletMapField
-  {
-  private:
-    Precision MinScalar;
-    Precision invDeltaScalar;
-
-  public:
-    VTKM_CONT
-    LerpScalar(const vtkm::Float32& minScalar, const vtkm::Float32& maxScalar)
-      : MinScalar(minScalar)
-    {
-      //Make sure the we don't divide by zero on
-      //something like an iso-surface
-      if (maxScalar - MinScalar != 0.f)
-        invDeltaScalar = 1.f / (maxScalar - MinScalar);
-      else
-        invDeltaScalar = 0.f;
-    }
-    typedef void ControlSignature(FieldIn<>,
-                                  FieldIn<>,
-                                  FieldIn<>,
-                                  FieldInOut<>,
-                                  WholeArrayIn<ScalarRenderingTypes>,
-                                  WholeArrayIn<>);
-    typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
-    template <typename ScalarPortalType, typename IndicesPortalType>
-    VTKM_EXEC void operator()(const vtkm::Id& hitIndex,
-                              const Precision& u,
-                              const Precision& v,
-                              Precision& lerpedScalar,
-                              const ScalarPortalType& scalars,
-                              const IndicesPortalType& indicesPortal) const
-    {
-      if (hitIndex < 0)
-        return;
-
-      vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
-
-      Precision n = 1.f - u - v;
-      Precision aScalar = Precision(scalars.Get(indices[1]));
-      Precision bScalar = Precision(scalars.Get(indices[2]));
-      Precision cScalar = Precision(scalars.Get(indices[3]));
-      lerpedScalar = aScalar * n + bScalar * u + cScalar * v;
-      //normalize
-      lerpedScalar = (lerpedScalar - MinScalar) * invDeltaScalar;
-    }
-  }; //class LerpScalar
-
-  template <typename Precision>
-  class NodalScalar : public vtkm::worklet::WorkletMapField
-  {
-  private:
-    Precision MinScalar;
-    Precision invDeltaScalar;
-
-  public:
-    VTKM_CONT
-    NodalScalar(const vtkm::Float32& minScalar, const vtkm::Float32& maxScalar)
-      : MinScalar(minScalar)
-    {
-      //Make sure the we don't divide by zero on
-      //something like an iso-surface
-      if (maxScalar - MinScalar != 0.f)
-        invDeltaScalar = 1.f / (maxScalar - MinScalar);
-      else
-        invDeltaScalar = 1.f / minScalar;
-    }
-
-    typedef void ControlSignature(FieldIn<>,
-                                  FieldOut<>,
-                                  WholeArrayIn<ScalarRenderingTypes>,
-                                  WholeArrayIn<>);
-
-    typedef void ExecutionSignature(_1, _2, _3, _4);
-    template <typename ScalarPortalType, typename IndicesPortalType>
-    VTKM_EXEC void operator()(const vtkm::Id& hitIndex,
-                              Precision& scalar,
-                              const ScalarPortalType& scalars,
-                              const IndicesPortalType& indicesPortal) const
-    {
-      if (hitIndex < 0)
-        return;
-
-      vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
-
-      //Todo: one normalization
-      scalar = Precision(scalars.Get(indices[0]));
-
-      //normalize
-      scalar = (scalar - MinScalar) * invDeltaScalar;
-    }
-  }; //class LerpScalar
-
-  template <typename Precision>
-  VTKM_CONT void Run(Ray<Precision>& rays,
-                     vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles,
-                     vtkm::cont::CoordinateSystem coordsHandle,
-                     const vtkm::cont::Field* scalarField,
-                     const vtkm::Range& scalarRange)
-  {
-    bool isSupportedField =
-      (scalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS ||
-       scalarField->GetAssociation() == vtkm::cont::Field::Association::CELL_SET);
-    if (!isSupportedField)
-      throw vtkm::cont::ErrorBadValue("Field not accociated with cell set or points");
-    bool isAssocPoints = scalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS;
-
-    // Find the triangle normal
-    vtkm::worklet::DispatcherMapField<CalculateNormals>(CalculateNormals())
-      .Invoke(
-        rays.HitIdx, rays.Dir, rays.NormalX, rays.NormalY, rays.NormalZ, coordsHandle, triangles);
-
-    // Calculate scalar value at intersection point
-    if (isAssocPoints)
-    {
-      vtkm::worklet::DispatcherMapField<LerpScalar<Precision>>(
-        LerpScalar<Precision>(vtkm::Float32(scalarRange.Min), vtkm::Float32(scalarRange.Max)))
-        .Invoke(rays.HitIdx, rays.U, rays.V, rays.Scalar, *scalarField, triangles);
-    }
-    else
-    {
-      vtkm::worklet::DispatcherMapField<NodalScalar<Precision>>(
-        NodalScalar<Precision>(vtkm::Float32(scalarRange.Min), vtkm::Float32(scalarRange.Max)))
-        .Invoke(rays.HitIdx, rays.Scalar, *scalarField, triangles);
-    }
-  } // Run
-
-}; // Class IntersectionData
-
-#define AABB_EPSILON 0.00001f
-class FindTriangleAABBs : public vtkm::worklet::WorkletMapField
-{
-public:
-  VTKM_CONT
-  FindTriangleAABBs() {}
-  typedef void ControlSignature(FieldIn<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                FieldOut<>,
-                                WholeArrayIn<Vec3RenderingTypes>);
-  typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
-  template <typename PointPortalType>
-  VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 4> indices,
-                            vtkm::Float32& xmin,
-                            vtkm::Float32& ymin,
-                            vtkm::Float32& zmin,
-                            vtkm::Float32& xmax,
-                            vtkm::Float32& ymax,
-                            vtkm::Float32& zmax,
-                            const PointPortalType& points) const
-  {
-    // cast to Float32
-    vtkm::Vec<vtkm::Float32, 3> point;
-    point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[1]));
-    xmin = point[0];
-    ymin = point[1];
-    zmin = point[2];
-    xmax = xmin;
-    ymax = ymin;
-    zmax = zmin;
-    point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[2]));
-    xmin = vtkm::Min(xmin, point[0]);
-    ymin = vtkm::Min(ymin, point[1]);
-    zmin = vtkm::Min(zmin, point[2]);
-    xmax = vtkm::Max(xmax, point[0]);
-    ymax = vtkm::Max(ymax, point[1]);
-    zmax = vtkm::Max(zmax, point[2]);
-    point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[3]));
-    xmin = vtkm::Min(xmin, point[0]);
-    ymin = vtkm::Min(ymin, point[1]);
-    zmin = vtkm::Min(zmin, point[2]);
-    xmax = vtkm::Max(xmax, point[0]);
-    ymax = vtkm::Max(ymax, point[1]);
-    zmax = vtkm::Max(zmax, point[2]);
-
-
-    vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
-    const vtkm::Float32 minEpsilon = 1e-6f;
-    xEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (xmax - xmin));
-    yEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (ymax - ymin));
-    zEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (zmax - zmin));
-
-    xmin -= xEpsilon;
-    ymin -= yEpsilon;
-    zmin -= zEpsilon;
-    xmax += xEpsilon;
-    ymax += yEpsilon;
-    zmax += zEpsilon;
-  }
-}; //class FindAABBs
-#undef AABB_EPSILON
-
-} // namespace detail
-
-class TriangleIntersector : public ShapeIntersector
+class VTKM_RENDERING_EXPORT TriangleIntersector : public ShapeIntersector
 {
 protected:
   vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Triangles;
   bool UseWaterTight;
 
 public:
-  TriangleIntersector()
-    : UseWaterTight(false)
-  {
-  }
+  TriangleIntersector();
 
-  void SetUseWaterTight(bool useIt) { UseWaterTight = useIt; }
+  void SetUseWaterTight(bool useIt);
 
   void SetData(const vtkm::cont::CoordinateSystem& coords,
-               vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles)
-  {
+               vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles);
 
-    CoordsHandle = coords;
-    Triangles = triangles;
-
-    vtkm::rendering::raytracing::AABBs AABB;
-    vtkm::worklet::DispatcherMapField<detail::FindTriangleAABBs>(detail::FindTriangleAABBs())
-      .Invoke(Triangles,
-              AABB.xmins,
-              AABB.ymins,
-              AABB.zmins,
-              AABB.xmaxs,
-              AABB.ymaxs,
-              AABB.zmaxs,
-              CoordsHandle);
-
-    this->SetAABBs(AABB);
-  }
-
-  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles() { return Triangles; }
-
-  class CellIndexFilter : public vtkm::worklet::WorkletMapField
-  {
-  public:
-    VTKM_CONT
-    CellIndexFilter() {}
-    typedef void ControlSignature(FieldInOut<>, WholeArrayIn<>);
-    typedef void ExecutionSignature(_1, _2);
-    template <typename TrianglePortalType>
-    VTKM_EXEC void operator()(vtkm::Id& hitIndex, TrianglePortalType& triangles) const
-    {
-      vtkm::Id cellIndex = -1;
-      if (hitIndex != -1)
-      {
-        cellIndex = triangles.Get(hitIndex)[0];
-      }
-
-      hitIndex = cellIndex;
-    }
-  }; //class CellIndexFilter
-
-  struct IntersectFunctor
-  {
-    template <typename Device, typename Precision>
-    VTKM_CONT bool operator()(Device,
-                              TriangleIntersector* self,
-                              Ray<Precision>& rays,
-                              bool returnCellIndex)
-    {
-      VTKM_IS_DEVICE_ADAPTER_TAG(Device);
-      self->IntersectRays(rays, Device(), returnCellIndex);
-      return true;
-    }
-  };
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles();
+  vtkm::Id GetNumberOfShapes() const override;
 
 
-  VTKM_CONT void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override
-  {
-    vtkm::cont::TryExecute(IntersectFunctor(), this, rays, returnCellIndex);
-  }
+  VTKM_CONT void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
+  VTKM_CONT void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
 
 
-  VTKM_CONT void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override
-  {
-    vtkm::cont::TryExecute(IntersectFunctor(), this, rays, returnCellIndex);
-  }
-
-  template <typename Precision, typename Device>
-  VTKM_CONT void IntersectRays(Ray<Precision>& rays,
-                               Device vtkmNotUsed(Device),
-                               bool returnCellIndex)
-  {
-
-    if (UseWaterTight)
-    {
-      WaterTightTriLeafIntersector<Device> leafIntersector(this->Triangles);
-
-      BVHTraverser<WaterTightTriLeafIntersector> traverser;
-      traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
-    }
-    else
-    {
-      MollerTriLeafIntersector<Device> leafIntersector(this->Triangles);
-
-      BVHTraverser<MollerTriLeafIntersector> traverser;
-      traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
-    }
-    // Normally we return the index of the triangle hit,
-    // but in some cases we are only interested in the cell
-    if (returnCellIndex)
-    {
-      vtkm::worklet::DispatcherMapField<CellIndexFilter> cellIndexFilterDispatcher;
-      cellIndexFilterDispatcher.Invoke(rays.HitIdx, Triangles);
-    }
-    // Update ray status
-    RayOperations::UpdateRayStatus(rays);
-  }
-
   VTKM_CONT void IntersectionData(Ray<vtkm::Float32>& rays,
                                   const vtkm::cont::Field* scalarField,
-                                  const vtkm::Range& scalarRange) override
-  {
-    IntersectionDataImp(rays, scalarField, scalarRange);
-  }
+                                  const vtkm::Range& scalarRange) override;
 
   VTKM_CONT void IntersectionData(Ray<vtkm::Float64>& rays,
                                   const vtkm::cont::Field* scalarField,
-                                  const vtkm::Range& scalarRange) override
-  {
-    IntersectionDataImp(rays, scalarField, scalarRange);
-  }
+                                  const vtkm::Range& scalarRange) override;
+
+  template <typename Precision>
+  VTKM_CONT void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
 
   template <typename Precision>
   VTKM_CONT void IntersectionDataImp(Ray<Precision>& rays,
                                      const vtkm::cont::Field* scalarField,
-                                     const vtkm::Range& scalarRange)
-  {
-    ShapeIntersector::IntersectionPoint(rays);
-    detail::TriangleIntersectionData intData;
-    intData.Run(rays, this->Triangles, this->CoordsHandle, scalarField, scalarRange);
-  }
+                                     const vtkm::Range& scalarRange);
 
-  vtkm::Id GetNumberOfShapes() const override { return Triangles.GetNumberOfValues(); }
 }; // class intersector
 }
 }