#include <petsc.h>

static PetscErrorCode ExtractStratumSubmesh(DM dm, PetscInt stratumValue)
{
  DM              subdm;
  DMLabel         FaceSetLabel, sublabel;
  IS              pointIS;
  const PetscInt *points;
  PetscBool       markedFaces = PETSC_TRUE;
  PetscInt        nitems;
  PetscInt        ibgn, iend;
  PetscMPIInt     rank, size;
  char            name[PETSC_MAX_PATH_LEN];

  PetscFunctionBeginUser;
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

  // Get the face set label
  PetscCall(DMGetLabel(dm, "Face Sets", &FaceSetLabel));
  PetscCheck(FaceSetLabel, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Face Sets returned a NULL label");

  PetscCall(DMLabelGetStratumIS(FaceSetLabel, stratumValue, &pointIS));
  PetscCheck(pointIS, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "DMLabelGetStratumIS returned a NULL label");
  PetscCall(ISGetSize(pointIS, &nitems));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "There are %" PetscInt_FMT " points in the Face Set with stratum %" PetscInt_FMT "\n", nitems, stratumValue));

  PetscCall(DMLabelCreate(PETSC_COMM_SELF, "boundary", &sublabel));
  PetscCall(ISGetIndices(pointIS, &points));
  for (PetscInt i = 0; i < nitems; ++i) PetscCall(DMLabelSetValue(sublabel, points[i], stratumValue));
  PetscCall(ISRestoreIndices(pointIS, &points));
  PetscCall(DMPlexLabelComplete(dm, sublabel));
  PetscCall(ISDestroy(&pointIS));

  PetscCall(DMPlexCreateSubmesh(dm, sublabel, stratumValue, markedFaces, &subdm));
  PetscCall(DMLabelDestroy(&sublabel));
  PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "patch_%2d.vtk", (int)stratumValue));
  PetscCall(PetscObjectSetName((PetscObject)subdm, name));
  PetscCall(DMPlexGetHeightStratum(subdm, 1, &ibgn, &iend));
  nitems = iend - ibgn;
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "There are %" PetscInt_FMT " point/faces in the Face Set with stratum %" PetscInt_FMT "\n", nitems, stratumValue));

  PetscCall(DMViewFromOptions(subdm, NULL, "-subdm_view"));
  PetscCall(DMDestroy(&subdm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char *argv[])
{
  DM dm;

  PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));

  PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
  PetscCall(DMSetType(dm, DMPLEX));
  PetscCall(DMSetFromOptions(dm));
  PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

  for (PetscInt i = 1; i <= 6; ++i) {
    PetscCall(ExtractStratumSubmesh(dm, i));
  }

  PetscCall(DMDestroy(&dm));
  PetscCall(PetscFinalize());
  return 0;
}