diff --git a/Packages/vcs/Lib/vcs2vtk.py b/Packages/vcs/Lib/vcs2vtk.py
index 31af7048de..88bd3adcbb 100644
--- a/Packages/vcs/Lib/vcs2vtk.py
+++ b/Packages/vcs/Lib/vcs2vtk.py
@@ -253,7 +253,7 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
 
     try:  # First try to see if we can get a mesh out of this
         g = data1.getGrid()
-        # Ok need unstrctured grid
+        # Ok need unstructured grid
         if isinstance(g, cdms2.gengrid.AbstractGenericGrid):
             continents = True
             wrap = [0., 360.]
@@ -263,7 +263,7 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
                 xM = m[:, 1].max()
                 ym = m[:, 0].min()
                 yM = m[:, 0].max()
-                N = m.shape[0]
+                numberOfCells = m.shape[0]
                 # For vtk we need to reorder things
                 m2 = numpy.ascontiguousarray(numpy.transpose(m, (0, 2, 1)))
                 m2.resize((m2.shape[0] * m2.shape[1], m2.shape[2]))
@@ -284,7 +284,8 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
                 xM = data2[:, 1].max()
                 ym = data2[:, 0].min()
                 yM = data2[:, 0].max()
-                N = data2.shape[0]
+                numberOfCells = data2.shape[0]
+                data2 = data2.filled(numpy.nan)
                 m2 = numpy.ascontiguousarray(numpy.transpose(data2, (0, 2, 1)))
                 nVertices = m2.shape[-2]
                 m2.resize((m2.shape[0] * m2.shape[1], m2.shape[2]))
@@ -299,20 +300,15 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
     if m3 is not None:
         # Create unstructured grid points
         vg = vtk.vtkUnstructuredGrid()
-        lst = vtk.vtkIdTypeArray()
-        cells = vtk.vtkCellArray()
-        numberOfCells = N
-        lst.SetNumberOfComponents(nVertices + 1)
-        lst.SetNumberOfTuples(numberOfCells)
-        for i in range(N):
-            tuple = [None] * (nVertices + 1)
-            tuple[0] = nVertices
+        for i in range(numberOfCells):
+            pt_ids = []
             for j in range(nVertices):
-                tuple[j + 1] = i * nVertices + j
-            lst.SetTuple(i, tuple)
-            # ??? TODO ??? when 3D use CUBE?
-        cells.SetCells(numberOfCells, lst)
-        vg.SetCells(vtk.VTK_POLYGON, cells)
+                indx = i * nVertices + j
+                if not numpy.isnan(m3[indx][0]):  # missing value means skip vertex
+                    pt_ids.append(indx)
+            vg.InsertNextCell(vtk.VTK_POLYGON,
+                              len(pt_ids),
+                              pt_ids)
     else:
         # Ok a simple structured grid is enough
         if grid is None:
diff --git a/testing/vcs/CMakeLists.txt b/testing/vcs/CMakeLists.txt
index 28e44a87d5..519c7dc655 100644
--- a/testing/vcs/CMakeLists.txt
+++ b/testing/vcs/CMakeLists.txt
@@ -23,6 +23,11 @@ cdat_add_test(vcs_test_png_set_size
   "${PYTHON_EXECUTABLE}"
   ${cdat_SOURCE_DIR}/testing/vcs/test_png_set_size.py
 )
+cdat_add_test(vcs_meshfill_missing_vertice
+  "${PYTHON_EXECUTABLE}"
+  ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_meshfill_vertices_issue.py
+  "${BASELINE_DIR}/test_vcs_meshfill_vertices_issue.png"
+)
 cdat_add_test(vcs_test_png_metadata
   "${PYTHON_EXECUTABLE}"
   ${cdat_SOURCE_DIR}/testing/vcs/test_png_metadata.py
@@ -153,7 +158,7 @@ cdat_add_test(vcs_test_background_mode_rotate
   ${cdat_SOURCE_DIR}/testing/vcs/test_background_mode_rotate.py
   )
 
-# Some vector plot testing
+# Some vector plot tests
 FOREACH(angle -180 -135 -90 -45 0 45 135 90 135 )
   cdat_add_test(vcs_test_vectors_angle_${angle}
     "${PYTHON_EXECUTABLE}"
diff --git a/testing/vcs/test_vcs_meshfill_vertices_issue.py b/testing/vcs/test_vcs_meshfill_vertices_issue.py
new file mode 100644
index 0000000000..6317ef0719
--- /dev/null
+++ b/testing/vcs/test_vcs_meshfill_vertices_issue.py
@@ -0,0 +1,57 @@
+import numpy
+import vcs
+import sys
+import os
+pth = os.path.join(os.path.dirname(__file__),"..")
+sys.path.append(pth)
+import checkimage
+
+x=vcs.init()
+
+x.setantialiasing(0)
+x.drawlogooff()
+x.setbgoutputdimensions(1200, 1090, units="pixels")
+
+data_values = [ 25, 45, 55.]
+
+data_lon = [ 5., 10., 15.]
+data_lat = [ 5., 10., 15.]
+
+data_lon_vert = [
+    # Triangle (last one missing because traingle has only 3 vertices
+    [2.5,7.5,5.,1.e20],
+    # Square
+    [7.5,12.5,12.5,7.5],
+    # Diamond
+    [15.,17.5,15,12.5],]
+
+data_lat_vert = [
+        # triangle
+        [2.5,2.5,7.5,1.e20],
+        # square
+        [7.5,7.5,12.5,12.5],
+        # diamond
+        [12.5,15,17.5,15],]
+
+mesh = numpy.array([data_lat_vert,data_lon_vert])
+print "MESH SHAPE:",mesh.shape,mesh.dtype
+mesh = numpy.transpose(mesh,(1,0,2))
+print "MESH SHAPE:",mesh.shape
+mesh = numpy.ma.masked_greater(mesh,1.e19)
+
+"""
+print "Triangle lats:"
+print mesh[0][0]
+print "Triangle lons:"
+print mesh[0][1]
+"""
+m = x.createmeshfill()
+m.levels = [20,30,50,70,80]
+m.mesh = True
+
+x.plot(numpy.array(data_values,),mesh,m,bg=True)
+x.png("test_vcs_meshfill_vertices_issue.png")
+src = sys.argv[1]
+ret = checkimage.check_result_image("test_vcs_meshfill_vertices_issue.png",
+                                    src, checkimage.defaultThreshold)
+sys.exit(ret)