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 <iostream>
#include <petsc.h>

//Connectivity of a 7x10 rectangular mesh of quads :
int 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):
int 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}
};

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  int            size, rank;
  ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

  if (size != 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a 2 processors example only");

  const int lLocalNumCells = 35; //Same on each rank for this example...
  const int lNumCorners = 4;
  const int lMeshDim    = 2;
  const int lGlobalNumVertices = 88;

  PetscInt lCells[lLocalNumCells*lNumCorners];

  for ( int i = 0; i < lLocalNumCells; ++i) {
    PetscInt* lNoeudCourant = lCells + i*lNumCorners;
    for ( int j = 0; j < 4; ++j ) {
      const int lElem = (sInitialPartition[rank])[i];
      (*lNoeudCourant) = sConnectivity[lElem][j];
      ++lNoeudCourant;
    }
  }

  DM       lDMBete = 0;

  ierr = DMPlexCreate(PETSC_COMM_WORLD,&lDMBete);CHKERRQ(ierr);

  ierr = DMSetDimension(lDMBete, lMeshDim);CHKERRQ(ierr);

  ierr = DMPlexBuildFromCellListParallel(lDMBete,
                                         lLocalNumCells,
                                         PETSC_DECIDE,
                                         lGlobalNumVertices,
                                         lNumCorners,
                                         lCells,
                                         PETSC_NULL); CHKERRQ(ierr);

  DM lDMBeteInterp = 0;

  ierr = DMPlexInterpolate(lDMBete, &lDMBeteInterp); CHKERRQ(ierr);
  ierr = DMDestroy(&lDMBete);  CHKERRQ(ierr);

  lDMBete = lDMBeteInterp;

  ierr = DMSetUseNatural(lDMBete,PETSC_TRUE);  CHKERRQ(ierr);

  const bool lDoThisToShowMeTheSegFault = true;
  PetscSection   section;
  if (lDoThisToShowMeTheSegFault) {
   PetscInt       numFields   = 1;
   PetscInt       numComp[1]  = {1};
   PetscInt       numDof[4]   = {0, 0, 1, 0}; //We create a unique field on elements to have the re-ordering done on a global vector containing element numbers....
   PetscInt       numBC       = 0;

   ierr = DMSetNumFields(lDMBete, numFields); CHKERRQ(ierr);

   ierr = DMPlexCreateSection(lDMBete, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section);  CHKERRQ(ierr);
   ierr = DMSetLocalSection(lDMBete, section); CHKERRQ(ierr);
  }

  PetscSF lSFMigrationSansOvl = 0;
  DM      lDMDistribueSansOvl = 0;
  PetscPartitioner lPart;
  DMPlexGetPartitioner(lDMBete, &lPart);
  PetscPartitionerSetFromOptions(lPart);

  ierr = DMPlexDistribute(lDMBete, 0, &lSFMigrationSansOvl, &lDMDistribueSansOvl); CHKERRQ(ierr);
  std::cout << "\n\n Migration SF:" << std::endl;

  ierr = PetscSFView(lSFMigrationSansOvl, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  PetscSF lSFPartition = 0;
  ierr = PetscSFCreateInverseSF(lSFMigrationSansOvl, &lSFPartition); CHKERRQ(ierr);

  std::cout << "\n\n Inversed Migration SF:" << std::endl;
  ierr = PetscSFView(lSFPartition, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

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

  PetscFinalize();
  return 0;
}