def field_subset(f, inds, rank=0):
    """Return the value of a field at a subset of points.

    Parameters
    ----------
    f: array, shape (a1, a2, ..., ad, r1, r2, ..., rrank)
        Rank-r field in d dimensions
    inds: integer array, shape (n, d)
        Index vectors
    rank: integer
        The rank of the field (0: scalar field, 1: vector field and so on).

    Returns
    -------
    f_sub: array, shape (n, rank)
        The subset of field values.
    """
    f_dim_space = f.ndim - rank
    if inds.ndim > 2:
        raise Exception('Too many dimensions in indices array')
    if inds.ndim == 1:
        if f_dim_space == 1:
            return f[inds]
        else:
            raise Exception('Indices array is 1d but field is not')
    if inds.shape[1] != f_dim_space:
        raise Exception('Indices and field dimensions do not match')
    # It's magic, don't touch it!
    return f[tuple([inds[:, i] for i in range(inds.shape[1])])]