program testDMView
  implicit none
#include 'finclude/petsc.h'
  DM :: dmMesh
  PetscErrorCode :: ierr
  MPI_Comm :: comm
  integer :: topologyDim
  integer :: numCell
  integer :: numNodePerCell
  integer :: geometryDim
  integer :: numNode
  integer,allocatable :: cellVertices(:)
  PetscReal, allocatable :: coordinates(:)
  PetscBool :: interpolate
  integer :: numX,numY
  PetscReal :: delX,delY
  integer :: ic,jc,kc
  integer :: numStrata
  IS :: cellIS ! cell id index set
  integer :: globalidStart,globalidEnd
  PetscReal :: v0(2),J(4),invJ(4),detJ
  integer :: cellId,coneSize
  PetscSection :: coordSection
  Vec :: CoordVec
  PetscScalar,pointer :: coords(:)
  integer :: locCrdSize
  ! Begin program
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
  comm = PETSC_COMM_SELF
  numX = 3
  numY = 3
  topologyDim=2
  numCell = (numX-1)*(numY-1)
  numNodePerCell = 4
  geometryDim = 2
  delX = 1.0/(numX-1)
  delY = 1.0/(numY-1)
  numNode = numX*numY
  allocate(cellVertices(numCell*numNodePerCell))
  allocate(coordinates(numNode*geometryDim))
  interpolate = .true.
  ! Prepare mesh Coordiantes
  kc=0
  do jc=1,numY
    do ic=1,numX
      kc=kc+1
      coordinates(kc) = (ic-1)*delX
      kc=kc+1
      coordinates(kc) = (jc-1)*delY
    end do
  end do
  ! Form Connectivity
  ! -1 becuase of zero based indexing for nodes
  kc=0
  do jc=1,numY-1
    do ic=1,numX-1
    kc=kc+1
    cellVertices((kc-1)*numNodePerCell + 1) = ic + (jc-1)*numX -1
    cellVertices((kc-1)*numNodePerCell + 2) = ic + (jc-1)*numX + 1 -1
    cellVertices((kc-1)*numNodePerCell + 3) = ic + (jc)*numX + 1 -1
    cellVertices((kc-1)*numNodePerCell + 4) = ic + (jc)*numX -1
    end do
  end do
  call DMPlexCreateFromCellList(comm,topologydim,numCell,numNode,numNodePerCell, &
                                interpolate,cellVertices,geometryDim,coordinates, &
                                dmMesh, ierr)
  print*,'Surcessfully created DMPlexSection'
  call DMView(dmMesh,PETSC_VIEWER_STDOUT_WORLD, ierr)

  call DMPlexGetStratumSize(dmMesh,"depth",1,numstrata,ierr)
  print*,'numstrate of vlaue',1,'is',numstrata
  call DMPlexGetStratumIS(dmMesh,"depth",1,cellIS,ierr)
  call ISView(cellIS,PETSC_VIEWER_STDOUT_WORLD, ierr)
  ! get the start and end ids of the 0-cells
  call DMPlexGetDepthStratum(dmMesh,0,globalidStart,globalidEnd,ierr)
  print*,'Node id start=',globalidStart,'node id end=',globalidEnd-1
  ! get the start and end ids of the 1-cells, 1-cell here is a edge
  call DMPlexGetDepthStratum(dmMesh,1,globalidStart,globalidEnd,ierr)
  print*,'edge id start=',globalidStart,'edge id end=',globalidEnd-1
  ! get the start and end ids of the 2-cells
  call DMPlexGetDepthStratum(dmMesh,2,globalidStart,globalidEnd,ierr)
  print*,'element id start=',globalidStart,'element id end=',globalidEnd-1
  ! compute cell geometry for first cell
  cellId = 0
  !-------------------------
  !<<<<<< No error here but Jacobian data seems weired >>>>>>>>>>>>>>>>>>>
  !------------------------
  call DMPlexComputeCellGeometry(dmMesh,cellId,v0,J,invJ,detJ,ierr)
  CHKERRQ(ierr)
  print*,'Translation Vector ::', v0
  print*,'Jacobian ::',J
  print*,'invJ ::',invJ
  print*,'detJ::',detJ
  print*,'Calling get ConeSize'
  !-------------------------
  !<<<<<< This call gives error >>>>>>>>>>>>>>>>>>>
  !------------------------
  call DMPlexGetConeSize(dmMesh, cellId,coneSize,ierr)
  CHKERRQ(ierr)
  call DMGetCoordinatesLocal(dmMesh,coordVec,ierr)
  !call VecView(coordVec,PETSC_VIEWER_STDOUT_WORLD, ierr)
  call DMPlexGetCoordinateSection(dmMesh, coordSection,ierr)
  !call PetscSectionView(coordSection, PETSC_VIEWER_STDOUT_WORLD, ierr)
  !-------------------------
  ! Now the problem begins, should not there be a DMPlexVecGetClosureF90 ?
  !<<<<<< This call gives error >>>>>>>>>>>>>>>>>>>
  !-------------------------
  call DMPlexVecGetClosure(dmMesh, coordSection, coordVec,cellId, loccrdSize,coords,ierr)
  !print*,coords


  call PetscFinalize(ierr)
end program testDMView