How to use the petsc4py.PETSc.IS function in petsc4py

def petsc_serial_matrix(test_space, trial_space, nnz=None):
    '''
    PETsc.Mat from trial_space to test_space to be filled in the 
    with block. The spaces can be represented by intergers meaning 
    generic R^n.
    '''
    # Decide local to global map
    # For our custom case everything is serial
    if is_number(test_space) and is_number(trial_space):
        comm = mpi_comm_world().tompi4py()
        # Local same as global
        sizes = [[test_space, test_space], [trial_space, trial_space]]

        row_map = PETSc.IS().createStride(test_space, 0, 1, comm)
        col_map = PETSc.IS().createStride(trial_space, 0, 1, comm)
    # With function space this can be extracted
    else:
        mesh = test_space.mesh()
        comm = mesh.mpi_comm().tompi4py()
        
        row_map = test_space.dofmap()
        col_map = trial_space.dofmap()
    
        sizes = [[row_map.index_map().size(IndexMap.MapSize_OWNED),
                  row_map.index_map().size(IndexMap.MapSize_GLOBAL)],
                 [col_map.index_map().size(IndexMap.MapSize_OWNED),
                  col_map.index_map().size(IndexMap.MapSize_GLOBAL)]]

        row_map = map(int, row_map.tabulate_local_to_global_dofs())
        col_map = map(int, col_map.tabulate_local_to_global_dofs())

self.pressureDOF = numpy.arange(start=L_range[0],
                                            stop=L_range[0]+neqns,
                                            step=nc,
                                            dtype="i")
            velocityDOF = []
            for start in range(1,nc):
                velocityDOF.append(numpy.arange(start=L_range[0]+start,
                                                stop=L_range[0]+neqns,
                                                step=nc,
                                                dtype="i"))
            self.velocityDOF = numpy.vstack(velocityDOF).transpose().flatten()
        self.pc = p4pyPETSc.PC().create()
        self.pc.setType('fieldsplit')
        self.isp = p4pyPETSc.IS()
        self.isp.createGeneral(self.pressureDOF,comm=p4pyPETSc.COMM_WORLD)
        self.isv = p4pyPETSc.IS()
        self.isv.createGeneral(self.velocityDOF,comm=p4pyPETSc.COMM_WORLD)
        self.pc.setFieldSplitIS(('velocity',self.isv),('pressure',self.isp))

rows.extend(shift + np.array(bc.keys(), dtype='int32'))
            x_values.extend(bc.values())
            
    rows = np.hstack(rows)
    x_values = np.array(x_values)

    x = bb.copy()
    x.zeroEntries()
    x.setValues(rows, x_values)

    # Apply to monolithic
    len(rows) and AA.zeroRowsColumns(rows, diag=diag_val, x=x, b=bb)

    blocks = []
    for first, last in zip(offsets[:-1], offsets[1:]):
        blocks.append(PETSc.IS().createStride(last-first, first, 1))

    # Reasamble
    comm = mpi_comm_world()
    b = [PETSc.Vec().createWithArray(bb.getValues(block), comm=comm) for block in blocks]
    for bi in b:
        bi.assemblyBegin()
        bi.assemblyEnd()
    b = block_vec(map(PETScVector, b))

    # PETSc 3.7.x
    try:
        AA.getSubMatrix
        A = [[PETScMatrix(AA.getSubMatrix(block_row, block_col)) for block_col in blocks]
             for block_row in blocks]
    # NOTE: 3.8+ does not ahve getSubMatrix, there is getLocalSubMatrix
    # but cannot get that to work - petsc error 73. So for now everything

dofmap = bc.function_space().dofmap()
            local_range = dofmap.ownership_range()
            local_dofs.update(list(range(local_range[0], local_range[1])))
            constrained_local_dofs.update([
                dofmap.local_to_global_index(local_dof_index) for local_dof_index in bc.get_boundary_values().keys()
            ])

        # List all unconstrained dofs
        unconstrained_local_dofs = local_dofs.difference(constrained_local_dofs)
        unconstrained_local_dofs = list(unconstrained_local_dofs)

        # Generate IS accordingly
        comm = bcs[0].function_space().mesh().mpi_comm()
        for bc in bcs:
            assert comm == bc.function_space().mesh().mpi_comm()
        self._is = PETSc.IS().createGeneral(unconstrained_local_dofs, comm)

def field_ises(self):
        """A list of PETSc ISes defining the global indices for each set in
        the DataSet.

        Used when extracting blocks from matrices for solvers."""
        ises = []
        nlocal_rows = 0
        for dset in self:
            nlocal_rows += dset.size * dset.cdim
        offset = self.comm.scan(nlocal_rows)
        offset -= nlocal_rows
        for dset in self:
            nrows = dset.size * dset.cdim
            iset = PETSc.IS().createStride(nrows, first=offset, step=1,
                                           comm=self.comm)
            iset.setBlockSize(dset.cdim)
            ises.append(iset)
            offset += nrows
        return tuple(ises)

step=3,
                                            dtype="i")
            velocityDOF = []
            for start in range(1,3):
                velocityDOF.append(numpy.arange(start=L_range[0]+start,
                                                stop=L_range[0]+neqns,
                                                step=3,
                                                dtype="i"))
                self.velocityDOF = numpy.vstack(velocityDOF).transpose().flatten()
        self.pc = p4pyPETSc.PC().create()
        if prefix:
            self.pc.setOptionsPrefix(prefix)
        self.pc.setType('fieldsplit')
        self.isp = p4pyPETSc.IS()
        self.isp.createGeneral(self.pressureDOF,comm=p4pyPETSc.COMM_WORLD)
        self.isv = p4pyPETSc.IS()
        self.isv.createGeneral(self.velocityDOF,comm=p4pyPETSc.COMM_WORLD)
        self.pc.setFieldSplitIS(('velocity',self.isv),('pressure',self.isp))
        self.pc.setFromOptions()
    def setUp(self, global_ksp=None):

if(i == i0n + In - 1 and j0 == j0n and k0 == k0n):
                ighosts += list(range(gindex+2*In*Jn+2*In*Kn+Jn*Kn,
                                 gindex+2*In*Jn+2*In*Kn+2*Jn*Kn))

        # zmax boundary
        i = i0+I
        if(i > Nz-1): i = 0
        for pd in pos_data:
            gindex, i0n, j0n, k0n, In, Jn, Kn = pd
            if(i == i0n and j0 == j0n and k0 == k0n):
                ighosts += list(range(gindex+2*In*Jn+2*In*Kn,
                                 gindex+2*In*Jn+2*In*Kn+Jn*Kn))

        # Set the ghost indices = finish constructing ghost vector
        self._ighosts = np.array(ighosts, dtype=np.int32)
        ghosts = PETSc.IS().createGeneral(self._ighosts)
        self._gvec.setMPIGhost(ghosts)