program main
      implicit none
#include "finclude/petsc.h90"
      
       
      PetscInt              :: N
      PetscErrorCode        :: ierr
      PetscInt, parameter   :: i1 = 1, i0 = 0
      PetscScalar,parameter :: r1 = 1.0, r0 = 0.0
      PetscBool             :: flg
      Vec                   :: globalVec
      Mat                   :: A
      DM                    :: dm
      integer               :: i,   t, myrank, comm_size, loc_vec_size,  volume
      PetscInt              :: pStart, pEnd, p, cStart, cEnd, c, vStart, vEnd, v, eStart, eEnd, e
      PetscInt              :: vecsize
      PetscSection          :: s
      real(8), dimension(:), pointer      :: ptr
    


      ! Setup PETSC
      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
      call MPI_Comm_rank(PETSC_COMM_WORLD,myrank,ierr)
      call MPI_Comm_size(PETSC_COMM_WORLD,comm_size,ierr)
      call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
      
      call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
      call DMPlexSetDimension(dm, 1*i1, ierr)

      !Consider following 1D mesh with periodic bounds. It consits of 3 boundaries (0,2), and
      !3 "volumes" (3,5)
      ! |--|--|--
      !
      ! So initalize 6 sieve points in total
      !
      call DMPlexSetChart(dm, i0, 6*i1, ierr)
      
      ! only the volumes have a cone which consits of the boundary points
      do volume = 3, 5 ! 
        call DMPlexSetConeSize(dm, volume*i1, 2*i1, ierr)
      enddo

      call DMSetUp(dm, ierr)

      ! tell plex the cone relations. periodic conditions show in last line

      call DMPlexSetCone(dm, 3*i1,[0*i1,1*i1],ierr)
      call DMPlexSetCone(dm, 4*i1,[1*i1,2*i1],ierr)
      call DMPlexSetCone(dm, 5*i1,[2*i1,0*i1],ierr)

      ! Now setup the supporting relationships for each sieve point
      ! That is reverse the arrows in the Hasse diagram and from which points
      ! p can now be reached (hope this description is correct)
      ! DMPLexSetSupportSize and DMPlexSetSupport exist, but instead of
      ! setting each point individually we just call DMPlexSymmetrize

      call DMPlexSymmetrize(dm,ierr)
      
      !In order to support efficient queries, we also want to construct fast 
      !search structures, indices, for the different
      !types of points, which is done using
      call DMPlexStratify(dm,ierr)

      ! This will give us the number of points in each level of the Hasse diag
      call DMPlexGetChart(dm, pStart, pEnd, ierr)
      call DMPlexGetHeightStratum(dm, 0, cStart, cEnd, ierr)
      call DMPlexGetHeightStratum(dm, 1, eStart, eEnd, ierr)
      call DMPlexGetDepthStratum(dm, 0, vStart, vEnd, ierr)

      call PetscSectionCreate(PETSC_COMM_WORLD, s, ierr)
      call PetscSectionSetChart(s, pStart, pEnd, ierr)

      print *, 'pstart', pstart, 'pEnd', pEnd 
      print *, 'cStart', cStart, 'cEnd', cEnd 
      print *, 'eStart', eStart, 'eEnd', eEnd
      print *, 'vStart', vStart, 'vEnd', vEnd


      ! Set the degrees of freedom
      ! Assume we want to calc. the transport of a scalar through the 1D mesh.
      ! So we have 1 dof on the "volumes"
      ! We also need transport coefficients on the boundaries. However, they
      ! should reside in a matrix => 0 dof on the boundaries
      ! 

      !do p = pStart, pEnd-1
        !call PetscSectionSetDof(s,p,1,ierr)
      !enddo
      do c = cStart, cEnd-1
        call PetscSectionSetDof(s,c,1,ierr)
      enddo
      !do v = vStart, vEnd-1
        !call PetscSectionSetDof(s,v,1,ierr)
      !enddo
      !do e = eStart, eEnd-1
        !call PetscSectionSetDof(s,e,0,ierr)
      !enddo

      call PetscSectionSetUp(s, ierr)

      ! Now lets get vectors!
      call DMSetDefaultSection(dm, s, ierr)
      call DMGetGlobalVector(dm, globalVec,ierr)

      call VecGetSize(globalVec,vecsize, ierr)
      print *, 'size of globalVec ', vecsize

      !Set the vector with contiguous numbers to see if we can retrieve
      !the right values using the plex routines

      do c = 0, vecsize-1 
        call VecSetValue(globalVec, c, c*r1, INSERT_VALUES, ierr)
      enddo
      call VecAssemblyBegin(globalVec,ierr)
      call VecAssemblyEnd(globalVec,ierr)
      call VecView(globalVec, PETSC_VIEWER_STDOUT_SELF, ierr)
  
      ! Why does is the closure of sieve point 4 not (1 , 2)?
      call DMPlexVecGetClosure(dm,s,globalVec,4*i1, ptr,ierr)
      print *, 'closure', ptr
     
      ! Choose FVM connection topology
      call DMPlexSetPreallocationCenterDimension(dm, i0, ierr);
      call DMCreateMatrix(dm, A, ierr)

      ! I guess know we would use DMPlexMatSetClosure to set the transport
      ! coefficients. However in the current setup the matrix contains 0
      ! writable elements :(

      call MatAssemblyBegin(A,ierr)
      call MatAssemblyEnd(A,ierr)

      call MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr)
      call MatView(A, PETSC_VIEWER_DRAW_WORLD, ierr)

      call DMDestroy(dm,ierr)

      call PETSCFINALIZE(ierr)
end program main