Fix Flying Edges z-boundary condition with non-square input

There was a bug with Flying Edges checking the boundary conditions in
the z-direction. It was comparing the z index to the size of the y
dimension, not the z dimension. This simple typo was probably missed
because most of the tests use square meshes. To catch the problem, you
need a mesh that has a different number of points in the y and z
directions and the contour has to go past the positive z boundary.

When this error was hit, the indices for the edges on that boundary were
not recorded, and later interpolation used garbage values for edge's
point ids.

data/baseline/filter/contour-uniform-boundaries.png

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:25d61b72591e1ecd66a70e1005cb31a5d37b0bba90bec42288dd41118801000e
+size 29986

docs/changelog/flying-edges-z-boundary-fix.md

@@ -0,0 +1,14 @@
+# Fixed Flying Edges Crash
+
+There was a bug in VTK-m's flying edges algorithm (in the contour filter
+for uniform structured data) that cause the code to read an index from
+uninitialized memory. This in turn caused memory reads from an
+inappropriate address that could cause bad values, failed assertions, or
+segmentation faults.
+
+The problem was caused by a misidentification of edges at the positive z
+boundary. Due to a typo, the z index was being compared to the length in
+the y dimension. Thus, the problem would only occur in the case where the y
+and z dimensions were of different sizes and the contour would go through
+the positive z boundary of the data, which was missing our test cases.
+

vtkm/filter/contour/testing/RenderTestContourFilter.cxx

@@ -10,11 +10,13 @@
 
 #include <vtkm/Math.h>
 #include <vtkm/cont/DataSet.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
 #include <vtkm/cont/testing/MakeTestDataSet.h>
 #include <vtkm/cont/testing/Testing.h>
 
 #include <vtkm/filter/clean_grid/CleanGrid.h>
 #include <vtkm/filter/contour/Contour.h>
+#include <vtkm/filter/field_transform/PointElevation.h>
 #include <vtkm/source/Tangle.h>
 
 #include <vtkm/io/VTKDataSetReader.h>
@@ -70,6 +72,40 @@ void TestContourFilterUniform()
     result, "pointvar", "filter/contour-uniform.png", testOptions);
 }
 
+void TestContourFilterUniformBoundaries()
+{
+  std::cout << "Generate Image for Contour filter on a uniform grid that goes through boundaries"
+            << std::endl;
+  // There was a bug in flying edges that did not identify boundaries when the dimension
+  // sizes were not the same.
+
+  vtkm::cont::DataSetBuilderUniform dsb;
+  vtkm::cont::DataSet dataSet =
+    dsb.Create({ 9, 5, 3 }, vtkm::Vec3f_64{ 0.0, 0.0, 0.0 }, vtkm::Vec3f_64{ 0.125, 0.25, 0.5 });
+
+  std::string fieldName = "pointvar";
+  vtkm::filter::field_transform::PointElevation elevation;
+  elevation.SetLowPoint(1.0, 0.0, 0.0);
+  elevation.SetHighPoint(0.0, 1.0, 1.0);
+  elevation.SetOutputFieldName(fieldName);
+  elevation.SetUseCoordinateSystemAsField(true);
+  dataSet = elevation.Execute(dataSet);
+
+  vtkm::filter::contour::Contour contour;
+  contour.SetGenerateNormals(true);
+  contour.SetMergeDuplicatePoints(true);
+  contour.SetIsoValues({ 0.25, 0.5, 0.75 });
+  contour.SetActiveField(fieldName);
+  vtkm::cont::DataSet result = contour.Execute(dataSet);
+
+  result.PrintSummary(std::cout);
+
+  //Y axis Flying Edge algorithm has subtle differences at a couple of boundaries
+  vtkm::rendering::testing::RenderTestOptions testOptions;
+  vtkm::rendering::testing::RenderTest(
+    result, fieldName, "filter/contour-uniform-boundaries.png", testOptions);
+}
+
 void TestContourFilterTangle()
 {
   std::cout << "Generate Image for Contour filter on a uniform tangle grid" << std::endl;
@@ -97,6 +133,7 @@ void TestContourFilterTangle()
 void TestContourFilter()
 {
   TestContourFilterUniform();
+  TestContourFilterUniformBoundaries();
   TestContourFilterTangle();
   TestContourFilterWedge();
 }

vtkm/filter/contour/worklet/contour/FlyingEdges.h

@@ -34,17 +34,7 @@ template <typename T, typename S>
 vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
 {
   vtkm::Id oldLen = handle.GetNumberOfValues();
-  if (oldLen == 0)
-  {
-    handle.Allocate(size);
-  }
-  else
-  {
-    vtkm::cont::ArrayHandle<T, S> tempHandle;
-    tempHandle.Allocate(oldLen + size);
-    vtkm::cont::Algorithm::CopySubRange(handle, 0, oldLen, tempHandle);
-    handle = tempHandle;
-  }
+  handle.Allocate(oldLen + size, vtkm::CopyFlag::On);
   return oldLen;
 }
 }

vtkm/filter/contour/worklet/contour/FlyingEdgesHelpers.h

@@ -15,6 +15,12 @@
 
 #include <vtkm/Types.h>
 
+#include <vtkm/cont/CellSetStructured.h>
+#include <vtkm/cont/DeviceAdapterTag.h>
+
+#include <vtkm/cont/cuda/internal/DeviceAdapterTagCuda.h>
+#include <vtkm/cont/kokkos/internal/DeviceAdapterTagKokkos.h>
+
 namespace vtkm
 {
 namespace worklet

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Common.h

@@ -185,7 +185,7 @@ struct Pass4TrimState
     {
       boundaryStatus[AxisToSum::zindex] += FlyingEdges3D::MinBoundary;
     }
-    if (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::yindex] - 2))
+    if (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2))
     {
       boundaryStatus[AxisToSum::zindex] += FlyingEdges3D::MaxBoundary;
     }

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h

@@ -15,8 +15,11 @@
 
 
 #include <vtkm/filter/contour/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Common.h>
 #include <vtkm/filter/contour/worklet/contour/FlyingEdgesTables.h>
 
+#include <vtkm/worklet/WorkletMapTopology.h>
+
 namespace vtkm
 {
 namespace worklet