def METHOD_NAME(ndim):
    """
    Test automatic averaging of Function at undefined grid points
    """
    grid = Grid(tuple([10]*ndim))
    dims = list(powerset(grid.dimensions))[1:]
    for d in dims:
        f = Function(name="f", grid=grid, staggered=d)
        # f at nod (x, y, z)
        shifted = f
        for dd in d:
            shifted = shifted.subs({dd: dd - dd.spacing/2})
        assert all(i == dd for i, dd in zip(shifted.indices, grid.dimensions))
        # Average automatically i.e.:
        # f not defined at x so f(x, y) = 0.5*f(x - h_x/2, y) + 0.5*f(x + h_x/2, y)
        avg = f
        for dd in d:
            avg = .5 * (avg + avg.subs({dd: dd - dd.spacing}))
        assert shifted.evaluate == avg