#include <petscdmplex.h>

PetscErrorCode SetVertexValues(DM dm, Vec *u)
{
  PetscScalar   *values;
  PetscInt       vStart, vEnd, v;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = DMCreateGlobalVector(dm, u);CHKERRQ(ierr);
  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
  ierr = VecGetArray(*u, &values);CHKERRQ(ierr);
  for (v = vStart; v < vEnd; ++v) {
    PetscScalar *p;

    ierr = DMPlexPointGlobalRef(dm, v, values, &p);CHKERRQ(ierr);
    if (p) p[0] = v + 1;
  }
  ierr = VecRestoreArray(*u, &values);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
  DM             dm, dmDist = NULL;
  Vec            v;
  PetscSection   section;
  IS             bcPointIS[1];
  PetscInt       dim = 2, numFields, numBC, i;
  PetscInt       numComp[3];
  PetscInt       numDof[12];
  PetscInt       bcField[1];
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);if (ierr) return ierr;
  ierr = PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);CHKERRQ(ierr);
  ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);
  ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr);
  if (dmDist) {
    ierr = DMDestroy(&dm);CHKERRQ(ierr);
    dm = dmDist;
  }
  numFields = 3;
  numComp[0] = 1;
  numComp[1] = dim;
  numComp[2] = dim - 1;
  for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0;
  numDof[0 * (dim + 1) + 0] = 1;
  numDof[1 * (dim + 1) + dim] = dim;
  numDof[2 * (dim + 1) + dim - 1] = dim - 1;
  numBC = 1;
  bcField[0] = 0;

  ierr = DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]);CHKERRQ(ierr);
  if (!bcPointIS[0]) numBC = 0;
  ierr = DMSetNumFields(dm, numFields);CHKERRQ(ierr);
  ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);CHKERRQ(ierr);
  ierr = ISDestroy(&bcPointIS[0]);CHKERRQ(ierr);
  /* Name the Field variables */
  ierr = PetscSectionSetFieldName(section, 0, "u");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 1, "v");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 2, "w");CHKERRQ(ierr);
  ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = DMSetSection(dm, section);CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);

  //DMGetPointSF(dm, &pointsf);
  //PetscSectionCreateGlobalSection(section, pointsf, PETSC_TRUE, PETSC_FALSE, &globalsection);
  //DMSetGlobalSection(dm, globalsection);
  ierr = SetVertexValues(dm, &v);CHKERRQ(ierr);
  ierr = VecView(v, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  //PetscSectionDestroy(&globalsection);
  ierr = VecDestroy(&v);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return (ierr);
}