#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;
    PetscViewer    viewer;
    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 = DMViewFromOptions(globalDM, NULL, "-dm_view");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 = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
    ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(viewer, "globalDM.vtu");CHKERRQ(ierr);
    ierr = DMView(globalDM, viewer);CHKERRQ(ierr);
    ierr = DMView(globalDM, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

    ierr = DMDestroy(&globalDM); CHKERRQ(ierr);

    ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
    ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(viewer, "subDM1.vtu");CHKERRQ(ierr);
    ierr = DMView(subDM[0], viewer);CHKERRQ(ierr);
    ierr = DMView(subDM[0], PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
    ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(viewer, "subDM2.vtu");CHKERRQ(ierr);
    ierr = DMView(subDM[1], viewer);CHKERRQ(ierr);
    ierr = DMView(subDM[1], PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
    ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(viewer, "ifacDM.vtu");CHKERRQ(ierr);
    ierr = DMView(ifacDM, viewer);CHKERRQ(ierr);
    ierr = DMView(ifacDM, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);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;
}