api: Use subs for origin in case index is a function

devito/operations/interpolators.py

@@ -189,8 +189,13 @@ def _interpolation_indices(self, variables, offset=0, field_offset=0,
                 mapper[d] = ConditionalDimension(p.name, self.sfunction._sparse_dim,
                                                  condition=condition, indirect=True)
 
+            # Apply mapper to each variable with origin correction before the
+            # Dimensions get replaced
+            subs = {v: v.subs({k: c - v.origin.get(k, 0) for k, c in mapper.items()})
+                    for v in variables}
+
             # Track Indexed substitutions
-            idx_subs.append(mapper)
+            idx_subs.append(subs)
 
         # Temporaries for the position
         temps = [Eq(v, k, implicit_dims=implicit_dims)

devito/tools/data_structures.py

@@ -66,6 +66,9 @@ def __getnewargs_ex__(self):
         # objects with varying number of attributes
         return (tuple(self), dict(self.__dict__))
 
+    def get(self, key, val):
+        return self._getters.get(key, val)
+
 
 class ReducerMap(MultiDict):
 

devito/types/basic.py

@@ -945,7 +945,8 @@ def origin(self):
         f(x) : origin = 0
         f(x + hx/2) : origin = hx/2
         """
-        return tuple(r - d for d, r in zip(self.dimensions, self.indices_ref))
+        return DimensionTuple(*(r-d for d, r in zip(self.dimensions, self.indices_ref)),
+                              getters=self.dimensions)
 
     @property
     def dimensions(self):
@@ -1197,8 +1198,8 @@ def indexify(self, indices=None, subs=None):
         subs = [{**{d.spacing: 1, -d.spacing: -1}, **subs} for d in self.dimensions]
 
         # Indices after substitutions
-        indices = [sympy.sympify((a - o).xreplace(s)) for a, o, s in
-                   zip(self.args, self.origin, subs)]
+        indices = [sympy.sympify(a.subs(d, d - o).xreplace(s)) for a, d, o, s in
+                   zip(self.args, self.dimensions, self.origin, subs)]
         indices = [i.xreplace({k: sympy.Integer(k) for k in i.atoms(sympy.Float)})
                    for i in indices]
         return self.indexed[indices]

devito/types/dense.py

@@ -1059,7 +1059,7 @@ def __indices_setup__(cls, **kwargs):
             mapper = {d: d for d in dimensions}
             for s in as_tuple(staggered):
                 c, s = s.as_coeff_Mul()
-                mapper.update({s: s + c * s.spacing/2})
+                mapper.update({s: s + c * s.spacing / 2})
             staggered_indices = mapper.values()
 
         return tuple(dimensions), tuple(staggered_indices)

tests/test_symbolics.py

@@ -9,7 +9,7 @@
                     Operator, SubDimension, norm, Le, Ge, Gt, Lt, Abs, sin, cos, Min, Max)
 from devito.ir import Expression, FindNodes
 from devito.symbolics import (retrieve_functions, retrieve_indexed, evalrel,  # noqa
-                              CallFromPointer, Cast, FieldFromPointer,
+                              CallFromPointer, Cast, FieldFromPointer, INT,
                               FieldFromComposite, IntDiv, ccode, uxreplace)
 from devito.types import Array, Bundle, LocalObject, Object, Symbol as dSymbol
 
@@ -31,6 +31,18 @@ def test_float_indices():
     assert indices == 1
 
 
+def test_func_of_indices():
+    """
+    Test that origin is correctly processed with functions
+    """
+    grid = Grid((10,))
+    x = grid.dimensions[0]
+    u = Function(name="u", grid=grid, space_order=2, staggered=x)
+    us = u.subs({u.indices[0]: INT(Abs(u.indices[0]))})
+    assert us.indices[0] == INT(Abs(x + x.spacing/2))
+    assert us.indexify().indices[0] == INT(Abs(x))
+
+
 @pytest.mark.parametrize('dtype,expected', [
     (np.float32, "float r0 = 1.0F/h_x;"),
     (np.float64, "double r0 = 1.0/h_x;")