examples: Boundary reconstruction

devito/operations/interpolators.py

@@ -80,13 +80,14 @@ class Injection(UnevaluatedSparseOperation):
     Evaluates to a list of Eq objects.
     """
 
-    def __new__(cls, field, expr, offset, interpolator, callback):
+    def __new__(cls, field, expr, offset, increment, interpolator, callback):
         obj = super().__new__(cls, interpolator, callback)
 
         # TODO: unused now, but will be necessary to compute the adjoint
         obj.field = field
         obj.expr = expr
         obj.offset = offset
+        obj.increment = increment
 
         return obj
 
@@ -255,7 +256,7 @@ def callback():
 
         return Interpolation(expr, offset, increment, self_subs, self, callback)
 
-    def inject(self, field, expr, offset=0):
+    def inject(self, field, expr, offset=0, increment=True):
         """
         Generate equations injecting an arbitrary expression into a field.
 
@@ -267,6 +268,8 @@ def inject(self, field, expr, offset=0):
             Injected expression.
         offset : int, optional
             Additional offset from the boundary.
+        increment: bool, optional
+            If True, generate increments (Inc) rather than assignments (Eq).            
         """
         def callback():
             # Derivatives must be evaluated before the introduction of indirect accesses
@@ -285,13 +288,18 @@ def callback():
                                                           field_offset=field_offset)
 
             # Substitute coordinate base symbols into the interpolation coefficients
-            eqns = [Inc(field.xreplace(vsub), _expr.xreplace(vsub) * b,
-                        implicit_dims=self.sfunction.dimensions)
-                    for b, vsub in zip(self._interpolation_coeffs, idx_subs)]
+            if increment:            
+                eqns = [Inc(field.xreplace(vsub), _expr.xreplace(vsub) * b,
+                            implicit_dims=self.sfunction.dimensions)
+                        for b, vsub in zip(self._interpolation_coeffs, idx_subs)]
+            else: 
+                eqns = [Eq(field.xreplace(idx_subs[0]), _expr.xreplace(idx_subs[0]),
+                implicit_dims=self.sfunction.dimensions)] 
+
 
             return temps + eqns
 
-        return Injection(field, expr, offset, self, callback)
+        return Injection(field, expr, offset, increment, self, callback)
 
 
 class PrecomputedInterpolator(GenericInterpolator):
@@ -388,3 +396,4 @@ def callback():
             return [Eq(_field, _field + rhs.subs(dim_subs))]
 
         return Injection(field, expr, offset, self, callback)
+

devito/types/sparse.py

@@ -848,7 +848,7 @@ def interpolate(self, expr, offset=0, u_t=None, p_t=None, increment=False):
                                                            increment=increment,
                                                            self_subs=subs)
 
-    def inject(self, field, expr, offset=0, u_t=None, p_t=None):
+    def inject(self, field, expr, offset=0, u_t=None, p_t=None, increment=True):
         """
         Generate equations injecting an arbitrary expression into a field.
 
@@ -864,14 +864,16 @@ def inject(self, field, expr, offset=0, u_t=None, p_t=None):
             Time index at which the interpolation is performed.
         p_t : expr-like, optional
             Time index at which the result of the interpolation is stored.
+        increment: bool, optional
+            If True, generate increments (Inc) rather than assignments (Eq).            
         """
         # Apply optional time symbol substitutions to field and expr
         if u_t is not None:
             field = field.subs({field.time_dim: u_t})
         if p_t is not None:
             expr = expr.subs({self.time_dim: p_t})
 
-        return super(SparseTimeFunction, self).inject(field, expr, offset=offset)
+        return super(SparseTimeFunction, self).inject(field, expr, offset=offset, increment=increment)
 
     # Pickling support
     _pickle_kwargs = AbstractSparseTimeFunction._pickle_kwargs +\

tests/test_interpolation.py

@@ -1,4 +1,4 @@
-from math import sin, floor
+from math import sin, floor, prod
 
 import numpy as np
 import pytest
@@ -53,6 +53,18 @@ def time_points(grid, ranges, npoints, name='points', nt=10):
     return points
 
 
+def time_grid_points(grid, name='points', nt=10):
+    """Create a SparseTimeFunction field with coordinates 
+       filled in by all grid points"""
+    npoints = prod(grid.shape) 
+    a = SparseTimeFunction(name=name, grid=grid, npoint=npoints, nt=nt)
+    dims = tuple([np.linspace(0., 1., d) for d in grid.shape]) 
+    for i in range(len(grid.shape)):
+        a.coordinates.data[:,i] = np.meshgrid(*dims)[i].flatten()
+
+    return a
+
+
 def a(shape=(11, 11)):
     grid = Grid(shape=shape)
     a = Function(name='a', grid=grid)
@@ -417,6 +429,31 @@ def test_inject_time_shift(shape, coords, result, npoints=19):
     assert np.allclose(a.data[indices], result, rtol=1.e-5)
 
 
+@pytest.mark.parametrize('shape, result, increment',  [
+    ((10, 10), 1., False),
+    ((10, 10), 5., True),
+    ((10, 10, 10), 1., False),
+    ((10, 10, 10), 5., True)        
+])
+def test_inject_time_increment(shape, result, increment):
+    """Test the increment option in the SparseTimeFunction's 
+       injection method. The increment=False option is 
+       expected to work only at points located on the grid, 
+       where no interpolation needed.
+    """
+    a = unit_box_time(shape=shape)
+    a.data[:] = 0.    
+    p = time_grid_points(a.grid, name='points', nt=10)
+
+    expr = p.inject(a, Float(1.), increment=increment)
+
+    Operator(expr)(a=a)    
+
+    assert np.allclose(a.data, result*np.ones(a.grid.shape), rtol=1.e-5)
+
+
+
+
 @pytest.mark.parametrize('shape, coords, result', [
     ((11, 11), [(.05, .95), (.45, .45)], 1.),
     ((11, 11, 11), [(.05, .95), (.45, .45), (.45, .45)], 0.5)