program main
#include <petsc/finclude/petscdmplex.h>
      use petscdmplex
      use petscsys
      implicit none

      DM :: dm, dmi
      PetscMPIInt       :: rank,ncpus
      PetscInt :: i
      PetscErrorCode ierr
      PetscInt  ::  numPoints(4), coneSize(20), cones(24)


!       9-------------8-------------13 
!      /|            /|            /|
!     / |           / |           / |
!    /  |          /  |          /  |
!   6-------------7-------------12  |
!   |   |         |   |         |   |
!   |   |         |   |         |   |
!   |   |         |   |         |   |
!   |   |         |   |         |   |
!   |   5------ --|---4---------|---11----------17----------16
!   |  /          |  /          |  /            /           /
!   | /           | /           | /            /           /
!   |/            |/            |/            /           /
!   2-------------3------------10-----------14----------15

      numPoints = (/ 0, 0, 0, 0 /)
      coneSize = (/ 8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4 /)
      cones = (/ 2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8,  10,14,17,11, 14,15,16,17 /)

      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
      if (ierr .ne. 0) then
        print*,'Unable to initialize PETSc'
        stop
      endif
      call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
      call MPI_Comm_size(PETSC_COMM_WORLD,ncpus,ierr)

      call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr)
      call PetscObjectSetName(dm, 'TestMesh', ierr);CHKERRA(ierr)
      call DMSetDimension(dm, 3, ierr);CHKERRA(ierr)
      call DMCreateLabel(dm, "celltype", ierr)

      if (rank == 0) then
        call DMPlexSetChart(dm, 0, 20, ierr);CHKERRA(ierr)
        do i=1,20
          call DMPlexSetConeSize(dm, i-1, coneSize(i), ierr)
        enddo
        call DMPlexSetCellType(dm, 0, DM_POLYTOPE_HEXAHEDRON, ierr)
        call DMPlexSetCellType(dm, 1, DM_POLYTOPE_HEXAHEDRON, ierr)
        call DMPlexSetCellType(dm, 18, DM_POLYTOPE_QUADRILATERAL, ierr)
        call DMPlexSetCellType(dm, 19, DM_POLYTOPE_QUADRILATERAL, ierr)
        do i=2,17
          call DMPlexSetCellType(dm, i, DM_POLYTOPE_POINT, ierr)
        enddo

        call DMSetUp(dm, ierr);CHKERRA(ierr)

        call DMPlexSetCone(dm, 0, cones(1:8), ierr);CHKERRA(ierr)
        call DMPlexSetCone(dm, 1, cones(9:16), ierr);CHKERRA(ierr)
        call DMPlexSetCone(dm, 18, cones(17:20), ierr);CHKERRA(ierr)
        call DMPlexSetCone(dm, 19, cones(21:24), ierr);CHKERRA(ierr)
      else
        call DMPlexSetChart(dm, 0, 0, ierr);CHKERRA(ierr)
      endif

      call DMPlexSymmetrize(dm,ierr);CHKERRA(ierr)
      call DMPlexStratify(dm,ierr);CHKERRA(ierr)
      call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)

      call DMPlexInterpolate(dm, dmi, ierr)
      call DMView(dmi, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)

      call DMDestroy(dm, ierr);CHKERRA(ierr)
      call DMDestroy(dmi, ierr);CHKERRA(ierr)
      call PetscFinalize(ierr)
    end program