gistfile1.py

import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc

bdy = PETSc.DMPlex().create()

bdy.setDimension(1)

bdy.createSquareBoundary([0, 0], [1, 1], [4, 4])

dm = PETSc.DMPlex().generate(bdy)

dm.markBoundaryFaces("boundary_faces")

dm.createLabel("ids")
coord_sec = dm.getCoordinateSection()
coords = dm.getCoordinates()
if dm.getStratumSize("boundary_faces", 1) > 0:
    for face in dm.getStratumIS("boundary_faces", 1).getIndices():
        c = dm.vecGetClosure(coord_sec, coords, face)
        if abs(c[0]) < 0.1 and abs(c[2]) < 0.1:
            dm.setLabelValue("ids", face, 1)
        if abs(c[0] - 1) < 0.1 and abs(c[2] - 1) < 0.1:
            dm.setLabelValue("ids", face, 2)
        if abs(c[1]) < 0.1 and abs(c[3]) < 0.1:
            dm.setLabelValue("ids", face, 3)
        if abs(c[1] - 1) < 0.1 and abs(c[3] - 1) < 0.1:
            dm.setLabelValue("ids", face, 4)

sf = dm.distribute()

dm.setRefinementUniform(True)

rdm = dm.refine()

# Grow overlap on refined mesh
sec = rdm.createSection([1], [1, 1, 1])
rdm.setDefaultSection(sec)
gsec = rdm.getDefaultGlobalSection()
sf = rdm.getDefaultSF()

lgmap = PETSc.LGMap().createSF(sf, gsec.getOffsetRange()[0])

rdm.distributeOverlap(lgmap, 1)

# Expect that we have all ids 1, ..., 4 covered (on at least one process)
for i in range(1, 5):
    if rdm.getStratumSize("ids", i) > 0:
        print PETSc.COMM_WORLD.rank, i, rdm.getStratumIS("ids", i).getIndices()