static char help[] = "Test FEM layout with constraints\n\n";

#include <petsc.h>

int main(int argc,char **argv) {
  DM                dm,pdm;
  PetscSection      section;
  const PetscInt    field = 0;
  char              ifilename[PETSC_MAX_PATH_LEN];
  PetscInt          sdim,s,pStart,pEnd,p,numVS,numPoints;
  PetscInt          constraints[1];
  IS                setIS,pointIS;
  const PetscInt   *setID,*pointID;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc,&argv,NULL,help));
  PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"FEM Layout Options","ex98");
  PetscCall(PetscOptionsString("-i","Filename to read","ex98",ifilename,ifilename,sizeof(ifilename),NULL));
  PetscOptionsEnd();

  PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD,ifilename,NULL,PETSC_TRUE,&dm));
  PetscCall(DMPlexDistributeSetDefault(dm,PETSC_FALSE));
  PetscCall(DMSetFromOptions(dm));

  PetscCall(DMPlexDistribute(dm,0,NULL,&pdm));
  if (pdm) {
    PetscCall(DMDestroy(&dm));
    dm = pdm;
  }
  PetscCall(DMViewFromOptions(dm,NULL,"-dm_view"));

  /* create a section */
  PetscCall(DMGetDimension(dm,&sdim));
  PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) dm),&section));
  PetscCall(PetscSectionSetNumFields(section,1));
  PetscCall(PetscSectionSetFieldName(section,field,"U"));
  PetscCall(PetscSectionSetFieldComponents(section,field,sdim));
  PetscCall(DMPlexGetChart(dm,&pStart,&pEnd));
  PetscCall(PetscSectionSetChart(section,pStart,pEnd));

  /* initialize the section storage for a P1 field */
  PetscCall(DMPlexGetDepthStratum(dm,0,&pStart,&pEnd));
  for (p = pStart; p < pEnd; ++p) {
    PetscCall(PetscSectionSetDof(section,p,sdim));
    PetscCall(PetscSectionSetFieldDof(section,p,0,sdim));
  }

  /* add constraints at all vertices belonging to a vertex set */
  /* first pass is to reserve space                            */
  PetscCall(DMGetLabelSize(dm,"Vertex Sets",&numVS));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,"# Vertex set: %" PetscInt_FMT "\n",numVS));
  PetscCall(DMGetLabelIdIS(dm,"Vertex Sets",&setIS));
  PetscCall(ISGetIndices(setIS,&setID));
  for (s = 0; s < numVS; ++s){
    PetscCall(DMGetStratumIS(dm,"Vertex Sets",setID[s],&pointIS));
    PetscCall(DMGetStratumSize(dm,"Vertex Sets",setID[s],&numPoints));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,"set %" PetscInt_FMT " size: %" PetscInt_FMT "\n",s,numPoints));
    PetscCall(ISGetIndices(pointIS,&pointID));
    for (p = 0; p < numPoints; ++p){
      PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\t point %" PetscInt_FMT "\n",pointID[p]));
      PetscCall(PetscSectionSetConstraintDof(section,pointID[p],1));
      PetscCall(PetscSectionSetFieldConstraintDof(section,pointID[p],field,1));
    }
    PetscCall(ISRestoreIndices(pointIS,&pointID));
    PetscCall(ISDestroy(&pointIS));
  }

  PetscCall(PetscSectionSetUp(section));

  /* add constraints at all vertices belonging to a vertex set          */
  /* second pass is to assign constraints to a specific component / dof */
  for (s = 0; s < numVS; ++s){
    PetscCall(DMGetStratumIS(dm,"Vertex Sets",setID[s],&pointIS));
    PetscCall(DMGetStratumSize(dm,"Vertex Sets",setID[s],&numPoints));
    PetscCall(ISGetIndices(pointIS,&pointID));
    for (p = 0; p < numPoints; ++p){
      constraints[0] = setID[s]%sdim;
      PetscCall(PetscSectionSetConstraintIndices(section,pointID[p],constraints));
      PetscCall(PetscSectionSetFieldConstraintIndices(section,pointID[p],field,constraints));
    }
    PetscCall(ISRestoreIndices(pointIS,&pointID));
    PetscCall(ISDestroy(&pointIS));
  }
  PetscCall(ISRestoreIndices(setIS,&setID));
  PetscCall(ISDestroy(&setIS));
  PetscCall(PetscObjectViewFromOptions((PetscObject) section,NULL,"-dm_section_view"));

  PetscCall(PetscSectionDestroy(&section));
  PetscCall(DMDestroy(&dm));

  PetscCall(PetscFinalize());
  return 0;
}

/*TEST
  build:
    requires: exodusii pnetcdf !complex
  testset:
    args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/SquareFaceSet.exo -dm_view -dm_section_view
    nsize: 1

    test:
      suffix: 0
      args:

TEST*/