Merge pull request #1691 from ccuetom/fix-precomputed-interpolation

API: fixed callback of PrecomputedInterpolator injection

devito/operations/interpolators.py

@@ -299,7 +299,7 @@ def __init__(self, obj, r, gridpoints_data, coefficients_data):
         gridpoints = SubFunction(name="%s_gridpoints" % self.obj.name, dtype=np.int32,
                                  dimensions=(self.obj.indices[-1], Dimension(name='d')),
                                  shape=(self._npoint, self.obj.grid.dim), space_order=0,
-                                 parent=self)
+                                 parent=self.obj)
 
         assert(gridpoints_data is not None)
         gridpoints.data[:] = gridpoints_data[:]
@@ -312,7 +312,7 @@ def __init__(self, obj, r, gridpoints_data, coefficients_data):
                                            shape=(self.obj.npoint, self.obj.grid.dim,
                                                   self.r),
                                            dtype=self.obj.dtype, space_order=0,
-                                           parent=self)
+                                           parent=self.obj)
         assert(coefficients_data is not None)
         interpolation_coeffs.data[:] = coefficients_data[:]
         self.obj._interpolation_coeffs = interpolation_coeffs
@@ -368,12 +368,12 @@ def callback():
             _expr = indexify(expr)
             _field = indexify(field)
 
-            p, _ = self.gridpoints.indices
+            p, _ = self.obj.gridpoints.indices
             dim_subs = []
             coeffs = []
-            for i, d in enumerate(self.grid.dimensions):
+            for i, d in enumerate(self.obj.grid.dimensions):
                 rd = DefaultDimension(name="r%s" % d.name, default_value=self.r)
-                dim_subs.append((d, INT(rd + self.gridpoints[p, i])))
+                dim_subs.append((d, INT(rd + self.obj.gridpoints[p, i])))
                 coeffs.append(self.obj.interpolation_coeffs[p, i, rd])
             rhs = prod(coeffs) * _expr
             _field = _field.subs(dim_subs)

tests/test_interpolation.py

@@ -169,6 +169,40 @@ def test_precomputed_interpolation_time():
         assert np.allclose(sf.data[it, :], it)
 
 
+def test_precomputed_injection():
+    """Test injection with PrecomputedSparseFunction which accepts
+       precomputed values for interpolation coefficients
+    """
+    shape = (11, 11)
+    coords = [(.05, .95), (.45, .45)]
+    origin = (0, 0)
+    result = 0.25
+
+    # Constant for linear interpolation
+    # because we interpolate across 2 neighbouring points in each dimension
+    r = 2
+
+    m = unit_box(shape=shape)
+    m.data[:] = 0.
+
+    gridpoints, interpolation_coeffs = precompute_linear_interpolation(coords,
+                                                                       m.grid, origin)
+
+    sf = PrecomputedSparseFunction(name='s', grid=m.grid, r=r, npoint=len(coords),
+                                   gridpoints=gridpoints,
+                                   interpolation_coeffs=interpolation_coeffs)
+
+    expr = sf.inject(m, FLOAT(1.))
+
+    Operator(expr)()
+
+    indices = [slice(0, 2, 1), slice(9, 11, 1)]
+    assert np.allclose(m.data[indices], result, rtol=1.e-5)
+
+    indices = [slice(4, 6, 1) for _ in coords]
+    assert np.allclose(m.data[indices], result, rtol=1.e-5)
+
+
 @pytest.mark.parametrize('shape, coords', [
     ((11, 11), [(.05, .9), (.01, .8)]),
     ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)])