#define PETSCDM_DLL
#include <petsc/private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/

#include <map>
#include <vector>
#include <set>

#include "petscdmplex.h"
#include "petscao.h"
#include "petscsnes.h"

int main(int argc, char *argv[])
{
    PetscErrorCode ierr;
    DM             globalDM;
    DM             distributedMesh = NULL;
    DM             subDM[2], ifacDM;
    PetscInt       dom1=-1, dom2=-1, ifac=-1;
    char           file[PETSC_MAX_PATH_LEN];
    PetscBool      set=PETSC_FALSE;
    DMLabel        cel_label, ifac_label;

    ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr);

    ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-Mesh", file, PETSC_MAX_PATH_LEN, &set);CHKERRQ(ierr);
    if(!set) {PetscFunctionReturn(1);}
    ierr = PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-domain1", &dom1, &set);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-domain2", &dom2, &set);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-interface", &ifac, &set);CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD, "Domain1: %d, Domain2: %d, Interface: %d.\n", dom1, dom2, ifac); CHKERRQ(ierr);
    
    ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, file, PETSC_TRUE, &globalDM);CHKERRQ(ierr);
    ierr = DMPlexDistribute(globalDM, 0, NULL, &distributedMesh);CHKERRQ(ierr);
    if (distributedMesh) {
        ierr = DMDestroy(&globalDM);CHKERRQ(ierr);
        globalDM = distributedMesh;
    }
    ierr = DMSetFromOptions(globalDM);CHKERRQ(ierr);

    ierr = DMGetLabel(globalDM, "Cell Sets", &cel_label); CHKERRQ(ierr);
    ierr = DMGetLabel(globalDM, "Face Sets", &ifac_label); CHKERRQ(ierr);
    ierr = DMLabelView(ifac_label, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

    ierr = DMPlexCreateSubmesh(globalDM, ifac_label, ifac, &ifacDM); CHKERRQ(ierr);
    ierr = DMPlexFilter(globalDM, cel_label, dom1, &subDM[0]); CHKERRQ(ierr);
    ierr = DMPlexFilter(globalDM, cel_label, dom2, &subDM[1]); CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD, "\n Global mesh \n"); CHKERRQ(ierr);
    ierr = DMViewFromOptions(globalDM, NULL, "-dm_view");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\n sub-mesh %d \n", dom1); CHKERRQ(ierr);
    ierr = DMViewFromOptions(subDM[0], NULL, "-dm_view");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\n sub-mesh %d \n", dom2); CHKERRQ(ierr);
    ierr = DMViewFromOptions(subDM[1], NULL, "-dm_view");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\n interface %d \n", ifac); CHKERRQ(ierr);
    ierr = DMViewFromOptions(ifacDM, NULL, "-dm_view");CHKERRQ(ierr);

    ierr = DMDestroy(&globalDM); CHKERRQ(ierr);
    ierr = DMDestroy(&subDM[0]); CHKERRQ(ierr);
    ierr = DMDestroy(&subDM[1]); CHKERRQ(ierr);
    ierr = DMDestroy(&ifacDM); CHKERRQ(ierr);
    
    ierr = PetscFinalize();CHKERRQ(ierr);
    
  return 0;
}