diff --git a/schism/geometry/skin.py b/schism/geometry/skin.py
index 5ce79c8..2e76a31 100644
--- a/schism/geometry/skin.py
+++ b/schism/geometry/skin.py
@@ -70,13 +70,13 @@ def _get_points(self):
 
         # Needs to use the interior mask of the target subgrid
         deriv_stagger = []
-        for dim in grid.dimensions:
+        for d in range(len(grid.dimensions)):
             try:
-                deriv_stagger.append(self.deriv.x0[dim])
+                deriv_stagger.append(self.deriv.x0[grid.dimensions[d]]
+                                     - grid.dimensions[d])
             except KeyError:
                 deriv_stagger.append(sp.core.numbers.Zero())
-        origin = tuple([self.deriv.expr.origin[d] - deriv_stagger[d]
-                        for d in range(len(grid.dimensions))])
+        origin = tuple(deriv_stagger)
 
         interior_mask = self.geometry.interior_mask[origin][all_points]
         mp = mp[:, interior_mask]
diff --git a/tests/test_skin.py b/tests/test_skin.py
index 54686a8..ead0bae 100644
--- a/tests/test_skin.py
+++ b/tests/test_skin.py
@@ -3,6 +3,7 @@
 import pytest
 import devito as dv
 import numpy as np
+import sympy as sp
 
 from schism.geometry.skin import stencil_footprint, ModifiedSkin
 
@@ -11,21 +12,29 @@ class DummyGeometry:
     """Dummy class to replace BoundaryGeometry for testing"""
     def __init__(self, grid):
         self.grid = grid
+        dims = grid.dimensions
         # Only works with square/cube grids
 
-        imask = np.zeros(grid.shape, dtype=bool)
+        imask = np.zeros((len(dims)+1,) + grid.shape, dtype=bool)
+        zero = sp.core.numbers.Zero()
 
         # Ugly but will do for now
         if len(self.grid.shape) == 2:  # 2D
+            origins = [(zero, zero), (dims[0].spacing/2, zero),
+                       (zero, dims[1].spacing/2)]
             bpoints = ([], [])
             for i in range(grid.shape[0]):
                 for j in range(grid.shape[1]):
                     if i == j:
                         bpoints[0].append(i)
                         bpoints[1].append(j)
+                        imask[1, i, j] = True
                     elif i > j:
-                        imask[i, j] = True
+                        imask[:, i, j] = True
         else:  # 3D
+            origins = [(zero, zero, zero), (dims[0].spacing/2, zero, zero),
+                       (zero, dims[1].spacing/2, zero),
+                       (zero, zero, dims[2].spacing/2)]
             bpoints = ([], [], [])
             for i in range(grid.shape[0]):
                 for j in range(grid.shape[1]):
@@ -34,13 +43,16 @@ def __init__(self, grid):
                             bpoints[0].append(i)
                             bpoints[1].append(j)
                             bpoints[2].append(k)
+                            imask[1, i, j, k] = True
                         elif i > k:
-                            imask[i, j, k] = True
+                            imask[:, i, j, k] = True
 
         self.boundary_points = tuple([np.array(bpoints[i])
                                       for i in range(len(bpoints))])
 
-        self.interior_mask = imask
+        self.interior_mask = {}
+        for i in range(len(origins)):
+            self.interior_mask[origins[i]] = imask[i]
 
 
 class TestMisc:
@@ -72,10 +84,15 @@ class TestSkin:
     """Tests for the ModifiedSkin object"""
     grid2d = dv.Grid(shape=(4, 4), extent=(10., 10.))
     grid3d = dv.Grid(shape=(4, 4, 4), extent=(10., 10., 10.))
+    x2d, y2d = grid2d.dimensions
+    x3d, y3d, z3d = grid3d.dimensions
 
     f2d = dv.TimeFunction(name='f', grid=grid2d, space_order=4)
     f3d = dv.TimeFunction(name='f', grid=grid3d, space_order=2)
 
+    v2d = dv.VectorTimeFunction(name='v', grid=grid2d, space_order=4)
+    v3d = dv.VectorTimeFunction(name='v', grid=grid3d, space_order=2)
+
     geom0 = DummyGeometry(grid2d)
     geom1 = DummyGeometry(grid3d)
 
@@ -87,12 +104,20 @@ class TestSkin:
     ans2 = (np.array([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]),
             np.array([0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]),
             np.array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]))
+    # Another answer
+    ans3 = (np.array([0, 1, 1, 2, 2, 2, 3, 3, 3]),
+            np.array([0, 0, 1, 0, 1, 2, 1, 2, 3]))
 
     @pytest.mark.parametrize('deriv, geometry, ans',
                              [(f2d.dx, geom0, ans0),
                               (f2d.dy, geom0, ans0),
                               (f2d.dxdy, geom0, ans1),
-                              (f3d.dx, geom1, ans2)])
+                              (f3d.dx, geom1, ans2),
+                              (v2d[0].dx(x0=x2d), geom0, ans0),
+                              (v2d[1].dy(x0=y2d), geom0, ans0),
+                              (v3d[0].dx(x0=x3d), geom1, ans2),
+                              (f2d.dx(x0=x2d+x2d.spacing/2), geom0, ans3),
+                              (f2d.dy(x0=y2d+y2d.spacing/2), geom0, ans0)])
     def test_modified_points(self, deriv, geometry, ans):
         """Check that modified points are correctly identified"""
         skin = ModifiedSkin(deriv, geometry)