program petscDMPlexTest
  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 :: interpolateFlag
  ! DMPlexsection related variables
  ! Becareful with the zero based indixening used
  ! in DMPlexsection input variables
  !integer :: tDim ! topologicalDim of the problem
  PetscSection ::   section
  integer :: numField ! number of Fields in the problem
  integer,allocatable :: numComp(:) ! size [numFiled]
  integer,allocatable :: numDof(:) ! size [numField*(tdim + 1)]
  ! the following is used to defines
  ! nodsets and facetsets for Boundary conditions
  integer :: numBC ! number of Boundary conditions
  integer,allocatable :: bcField(:) ! size [numBc]
  IS,allocatable :: bcPointIS(:) ! size [numBC], each IS holds sieve points to which BC applies
  integer,allocatable :: numDofperField(:,:) ! size [tdim+1,numField]
  integer :: tdim, FieldId ! loop variables

  ! Begin program
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)

  comm =  PETSC_COMM_SELF
  topologyDim = 2 !
  interpolateFlag = .true.
  ! This line gives error. Looks like there is no Fortran interface
  call DMPlexCreateBoxMesh(comm, topologyDim, interpolateFlag, dmMesh, ierr)
  !-----------------------------------------------------
  ! Convention for topology dimension
  ! point or node has topology dim 0
  ! line or edge has topology dim 1
  ! Traingle or qudrangle has topology dim 2
  ! Tetrahedral or other three dimensional element has topology dim 3
  !------------------------------------------------------
  ! Temporary definiton of varibles, change this to single scalr variable 
  ! later to test laplace problem
  ! Create a section assuming one scalar fields u and one vector field v
  numField = 2
  allocate(numComp(numField))
  numComp = [1,2] ! since u is scalar and v is two component vector (mesh has geometric dimension 2)
  ! Square mesh so tdim = 2
  ! Let the scalar field u be deinfed on nodes and v on cell
  allocate(numDofperField(0:topologyDim,0:numField-1))
  numDofperField(:,0) = [4,0,0] ! Field variable u : 4 is the number of nodes per cell, quadrangle elements
  numDofperField(:,1) = [0,0,1] ! Field Vraiable v : 1 is the number of cell
  allocate(numDof(numField*(topologyDim+1)))
  ! collect all Dofs in the cell
  ! basically just reshaping numDofperField to a 1D array
  do fieldId=0,numField-1,1 ! zero based indexing
    do tdim=0,topologyDim,1
      numDof(fieldId*(topologyDim+1)+tdim) = numDofperField(tdim,fieldId)
    end do
  end do
  numBC = 2 ! Two Dirchilet Boundary conditons
  allocate(bcField(numBC))
  bcField = [0,0] ! both boundary conditions on field u
  allocate(bcPointIS(numBC))
  ! get the index set for each marker
  ! I do not know what stratum value is, used 1 here
  ! Ask Mathew for help
  ! Idea is to set Dirichelet BC on left an right boundary of 
  ! a suqare box mesh: (0,0) and (1,1) 
  ! left boundary is x=0, y = 0 to 1.0
  ! right boundary is x=1.0, y = 0 to 1.0
  call DMPlexGetStratumIS(dmMesh, "source", 1, bcPointIS(1),ierr)
  call DMPlexGetStratumIS(dmMesh, "drain", 1, bcPointIS(2),ierr)
  !
  call DMPlexCreateSection(dmMesh, topologyDim, numField, numComp, numDof, &
                           numBC, bcField, bcPointIS,section,ierr)
  ! Name the Field variables
  call PetscSectionSetFieldName(section, 0, "scalarNodeVar",ierr) ! 0 is field id
  call PetscSectionSetFieldName(section, 1, "vectorCellVar",ierr) ! 1 is field id
  ! Copied form ex62.c no idea what this is
  call DMSetDefaultSection(dmMesh, section,ierr)
  print*,'Surcessfully created DMPlexSection'
  call PetscFinalize(ierr)
end program petscDMPlexTest