How to use the pytools.indices_in_shape function in pytools

def split(self, whole_vol_vector):
        from pytools import indices_in_shape
        from hedge.tools import log_shape

        ls = log_shape(whole_vol_vector)

        if ls != ():
            result = [numpy.zeros(ls, dtype=object)
                    for part_emb in self._embeddings()]
            for p, part_emb in enumerate(self._embeddings()):
                for i in indices_in_shape(ls):
                    result[p][i] = whole_vol_vector[part_emb]
            return result
        else:
            return [whole_vol_vector[part_emb]
                    for part_emb in self._embeddings()]

def apply_index_map(imap, vector):
    from hedge._internal import VectorTarget, perform_index_map

    ls = log_shape(vector)
    if ls == ():
        result = numpy.zeros(shape=(imap.to_length,), dtype=float)
        perform_index_map(imap, VectorTarget(vector, result))
    else:
        result = numpy.zeros(shape=ls+(imap.to_length,), dtype=float)
        from pytools import indices_in_shape
        for i in indices_in_shape(ls):
            perform_index_map(imap, VectorTarget(vector[i], result[i]))
    return result

def with_object_array_or_scalar(f, field):
    ls = log_shape(field)
    if ls != ():
        from pytools import indices_in_shape
        result = numpy.zeros(ls, dtype=object)
        for i in indices_in_shape(ls):
            result[i] = f(field[i])
        return result
    else:
        return f(field)

"""

        from pytools import indices_in_shape
        from scipy.special import eval_chebyt

        def eval_basis(ind, x):
            result = 1
            for i in range(self.dim):
                coord = (x[i] - self.center[i])/(self.h/2)
                result *= eval_chebyt(ind[i], coord)
            return result

        from functools import partial
        return [
                partial(eval_basis, ind)
                for ind in indices_in_shape((self.npoints,)*self.dim)]

def volume_from_gpu(self, field, check=None):
        if check is None:
            check = self.debug

        from hedge.tools import log_shape
        ls = log_shape(field)
        if ls != ():
            result = numpy.array(ls, dtype=object)

            from pytools import indices_in_shape

            for i in indices_in_shape(ls):
                result[i] = self.from_gpu(field[i])
            return result
        else:
            copied_vec = field.get(pagelocked=True)
            result = numpy.empty(shape=(len(self),), dtype=copied_vec.dtype)

            block_dofs = self.block_dof_count()

            block_offset = 0
            for block in self.blocks:
                face_length = block.local_discretization.face_node_count()
                el_length = block.local_discretization.node_count()

                # write internal dofs
                el_offset = block_offset
                for cpu_slice in block.cpu_slices:

def _new_vec(self, shape, create_func, dtype, base_size):
        if dtype is None:
            dtype = self.default_scalar_type

        if shape == ():
            return create_func((base_size,), dtype=dtype)

        result = numpy.empty(shape, dtype=object)
        from pytools import indices_in_shape
        for i in indices_in_shape(shape):
            result[i] = create_func((base_size,), dtype=dtype)
        return result

def reassemble(self, parts_vol_vectors):
        from pytools import single_valued, indices_in_shape
        from hedge.tools import log_shape
        ls = single_valued(log_shape(pvv) for pvv in parts_vol_vectors)

        def remap_scalar_field(idx):
            result = self.whole_discr.volume_zeros()
            for part_emb, part_vol_vector in zip(
                    self._embeddings(), parts_vol_vectors):
                result[part_emb] = part_vol_vector[idx]

            return result

        if ls != ():
            result = numpy.zeros(ls, dtype=object)
            for i in indices_in_shape(ls):
                result[i] = remap_scalar_field(i)
            return result
        else:
            return remap_scalar_field(())

def __call__(self, from_vec):
        from hedge._internal import perform_elwise_operator
        from hedge.tools import log_shape

        ls = log_shape(from_vec)
        result = np.empty(shape=ls, dtype=object)

        from pytools import indices_in_shape
        for i in indices_in_shape(ls):
            result_i = self.to_discr.volume_zeros(kind="numpy")
            result[i] = result_i

            for from_eg, to_eg, imat in zip(
                    self.from_discr.element_groups,
                    self.to_discr.element_groups,
                    self.interp_matrices):
                perform_elwise_operator(
                        from_eg.ranges, to_eg.ranges,
                        imat, from_vec[i], result_i)

        if ls == ():
            return result[()]
        else:
            return result

def _new_bdry(self, tag, shape, create_func):
        result = numpy.empty(shape, dtype=object)
        from pytools import indices_in_shape
        bdry = self.get_boundary(TAG_ALL)
        for i in indices_in_shape(shape):
            result[i] = create_func(
                    (self.aligned_boundary_floats), 
                    dtype=numpy.float32)
        return result