static char help[] = "The main goal of this code is to retrieve the original element numbers as found in the "
                     "initial partitions (sInitialPartition)... but after the call to DMPlexDistribute";

#include <petsc.h>

//Connectivity of a 7x10 rectangular mesh of quads :
PetscInt sConnectivity[70][4] = {
{0,4,34,22},
{22,34,35,23},
{23,35,36,24},
{24,36,37,25},
{25,37,38,26},
{26,38,39,27},
{27,39,13,2},
{4,5,40,34},
{34,40,41,35},
{35,41,42,36},
{36,42,43,37},
{37,43,44,38},
{38,44,45,39},
{39,45,14,13},
{5,6,46,40},
{40,46,47,41},
{41,47,48,42},
{42,48,49,43},
{43,49,50,44},
{44,50,51,45},
{45,51,15,14},
{6,7,52,46},
{46,52,53,47},
{47,53,54,48},
{48,54,55,49},
{49,55,56,50},
{50,56,57,51},
{51,57,16,15},
{7,8,58,52},
{52,58,59,53},
{53,59,60,54},
{54,60,61,55},
{55,61,62,56},
{56,62,63,57},
{57,63,17,16},
{8,9,64,58},
{58,64,65,59},
{59,65,66,60},
{60,66,67,61},
{61,67,68,62},
{62,68,69,63},
{63,69,18,17},
{9,10,70,64},
{64,70,71,65},
{65,71,72,66},
{66,72,73,67},
{67,73,74,68},
{68,74,75,69},
{69,75,19,18},
{10,11,76,70},
{70,76,77,71},
{71,77,78,72},
{72,78,79,73},
{73,79,80,74},
{74,80,81,75},
{75,81,20,19},
{11,12,82,76},
{76,82,83,77},
{77,83,84,78},
{78,84,85,79},
{79,85,86,80},
{80,86,87,81},
{81,87,21,20},
{12,1,28,82},
{82,28,29,83},
{83,29,30,84},
{84,30,31,85},
{85,31,32,86},
{86,32,33,87},
{87,33,3,21}
};

//The initial partitions given by reading (simulating a read by blocks for large meshes):
PetscInt sInitialPartition[2][35] = {
  {0,1,2,6,7,8,12,13,14,18,19,20,24,25,26,30,31,32,36,37,38,42,43,44,48,49,50,54,55,56,60,61,62,66,67},
  {3,4,5,9,10,11,15,16,17,21,22,23,27,28,29,33,34,35,39,40,41,45,46,47,51,52,53,57,58,59,63,64,65,68,69}
};

int main(int argc, char **argv)
{
  const PetscInt   Nc   = 35; //Same on each rank for this example...
  const PetscInt   Ncor = 4;
  const PetscInt   dim  = 2;
  const PetscInt   Nv   = 88;
  DM               dm, idm, ddm;
  PetscSF          sfVert, sfMig, sfPart;
  PetscPartitioner part;
  PetscSection     s;
  PetscInt        *cells, c;
  PetscMPIInt      size, rank;
  PetscBool        box = PETSC_FALSE, field = PETSC_FALSE;
  PetscErrorCode   ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr);
  if (size != 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a 2 processors example only");
  ierr = PetscOptionsGetBool(NULL, NULL, "-box", &box, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL, NULL, "-field", &field, NULL);CHKERRQ(ierr);

  ierr = DMPlexCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr);
  if (box) {
    ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
  } else {
    ierr = PetscMalloc1(Nc * Ncor, &cells);CHKERRQ(ierr);
    for (c = 0; c < Nc; ++c) {
      PetscInt cell = sInitialPartition[rank][c], cor;

      for (cor = 0; cor < Ncor; ++cor) {
        cells[c*Ncor + cor] = sConnectivity[cell][cor];
      }
    }
    ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
    ierr = DMPlexBuildFromCellListParallel(dm, Nc, PETSC_DECIDE, Nv, Ncor, cells, &sfVert);CHKERRQ(ierr);
    //ierr = DMPlexBuildCoordinatesFromCellListParallel(dm, dim, sfVert, coords);CHKERRQ(ierr);
    ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);
    ierr = PetscFree(cells);CHKERRQ(ierr);
    ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
    ierr = DMDestroy(&dm);CHKERRQ(ierr);
    dm   = idm;
  }
  ierr = DMSetUseNatural(dm, PETSC_TRUE);CHKERRQ(ierr);
  ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);

  if (field) {
   const PetscInt Nf         = 1;
   const PetscInt numComp[1] = {1};
   const PetscInt numDof[3]  = {1, 0, 0};
   const PetscInt numBC      = 0;

   ierr = DMSetNumFields(dm, Nf); CHKERRQ(ierr);
   ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &s);  CHKERRQ(ierr);
   ierr = DMSetLocalSection(dm, s); CHKERRQ(ierr);
   ierr = PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
  }

  ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr);
  ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);

  ierr = DMPlexDistribute(dm, 0, &sfMig, &ddm); CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscSFView(sfMig, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  ierr = PetscSFCreateInverseSF(sfMig, &sfPart); CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) sfPart, "Inverse Migration SF");CHKERRQ(ierr);
  ierr = PetscSFView(sfPart, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  //... to be continued....

  ierr = PetscSFDestroy(&sfMig);CHKERRQ(ierr);
  ierr = PetscSFDestroy(&sfPart);CHKERRQ(ierr);
  ierr = DMDestroy(&ddm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

/*TEST

  testset:
    args: -field

    test:
      suffix: 0
      args:

    test:
      suffix: 1
      args: -box -dm_plex_simplex 0 -dm_plex_box_faces 7,10 -dm_distribute -dm_view hdf5:mesh.h5

TEST*/