program testDMView
  implicit none
#include 'finclude/petsc.h90'
  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,dim
  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
  integer :: facetId
  integer,pointer :: sharedCellId(:)
  ! Begin program
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
  comm = PETSC_COMM_SELF
  !------------------------------
  ! Create a box mesh with quadrangle elements
  ! Mesh region = [0 1] x [0 1]
  ! numX -  number of points in x-direction
  ! numY -  number of points in y-direction
  !------------------------------
  numX = 3
  numY = 3
  ! topology dimension is 2 since mesh is two dimentsional
  topologyDim=2
  ! Number of Cells
  numCell = (numX-1)*(numY-1)
  ! Number of Nodes per Cell is 4 for quadrangle element
  numNodePerCell = 4
  ! geometry dimensio is 2, since the mesh is in x-y plane
  geometryDim = 2
  ! Grid Spacing in x and y direction
  delX = 1.0/(numX-1)
  delY = 1.0/(numY-1)
  ! Number of Nodes/vertices in the mesh
  numNode = numX*numY
  allocate(cellVertices(numCell*numNodePerCell))
  allocate(coordinates(numNode*geometryDim))
  interpolate = .true.
  ! Prepare mesh Coordiantes
  ! Mesh Nodes are numberd in x-direction first then in y direction
  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
  ! Cells are counted in X-direction first then in Y-direction
  ! In each cells the nodes are listed in counter clockwise direction
  ! For example, for cellid = 0, the node coordiantes are
  ! node1 = [0,0] node2 = [0,delX], node3 = [delX,delY], node4 = [0,delY]
  !kc=0
  do jc=1,numY-1
    do ic=1,numX-1
    kc=ic + (numX-1)*(jc-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
  !--------------------------------
  ! Create DM using vertices and Connectivity
  ! Any entity of a mesh is identified as cell with dimension dim
  ! or a cell with codimension cdim ?
  ! Terminology Notes:
  ! vertex = 0 dimensional cell or dim codimension cell
  ! element = dim dimensional cell or 0 codimension cell
  ! facet = dim-1 dimensional cell
  !
  ! node or point = 0 dimension topology
  ! edge or line  = 1 dimension topology
  ! Triangle or qudrangle  = 2 dimension topology
  ! Terahedron or pyramid or Any 3d element is = 3 dimension topology
  ! By defualt DM created from cell list has only 0 and dim dimensional entities
  ! use interpolateflag=.true. to create other entites
  !--------------------------------
  call DMPlexCreateFromCellList(comm,topologydim,numCell,numNode,numNodePerCell, &
                                interpolate,cellVertices,geometryDim,coordinates, &
                                dmMesh, ierr)
  print*,'Surcessfully created DM From Cell List'
  !--------------------------------
  ! Lets print the dmMesh using DMView to see how the DM information
  ! is organized
  !--------------------------------
  call DMView(dmMesh,PETSC_VIEWER_STDOUT_WORLD, ierr)
  !-------------------------------
  ! You should see an output like this on screen
  ! Mesh in 2 dimensions:
  !   0-cells: 9  ! = numX*numY
  !   1-cells: 12
  !   2-cells: 4  ! = (numX-1)*(numY-1)
  ! Labels:
  !   depth: 3 strata of sizes (9, 12, 4)
  !--------------------------------

  ! get the start and end ids of the 1-cells, 1-cell here is a edge
  print*,'printing support for facet'
  call DMPlexGetDepthStratum(dmMesh,1,globalidStart,globalidEnd,ierr)
  print*,'edge id start=',globalidStart,'edge id end=',globalidEnd-1
  do facetId=globalidStart,globalidEnd-1
  call DMPlexGetSupport(dmMesh,facetId,SharedCellId,ierr)
  print*,'f:',facetId,'c:',sharedCellId
  call DMPlexRestoreSupport(dmMesh,facetId,SharedCellId,ierr)
  end do
  ! Now print the support for each facet

  call PetscFinalize(ierr)
end program testDMView