Was there any particular reason for making DMPlexMarkBoundaryFaces a quadratic time algorithm? I realize it's hard to write a subquadratic time algorithm on top of DMLabelSetValue; maybe PETSc needs some basic integer hash tables?
Note: boundaries are not necessarily sublinear in the case of strongly adapted meshes. Geoffrey
