program test

#include "petsc/finclude/petscvec.h"
#include "petsc/finclude/petscdmplex.h"
#include "petsc/finclude/petscdmlabel.h"
 use petscvec
 use petscdmplex
 use petscsys

 implicit none

 DM                                :: dm_dumm, dm_geom
 Vec                               :: vglo, vnat, coord
 PetscSF                           :: sf_nat
 PetscMPIInt                       :: myid, nproc
 PetscErrorCode                    :: ierr
 PetscBool                         :: Check
 PetscInt                          :: dim
 DMLabel, pointer                  :: nolabel(:) => NULL()
 PetscSection                      :: section
 PetscInt                          :: numFields = 3 
 PetscInt                          :: numBC     = 1
 PetscInt, target, dimension(3)    :: numComp 
 PetscInt, target, dimension(12)   :: numDof  
 PetscInt, pointer                 :: pNumComp(:)
 PetscInt, pointer                 :: pNumDof(:)
 PetscInt, target, dimension(1)    :: bcField
 PetscInt, pointer                 :: pBcField(:)
 IS, target, dimension(1)          :: bcCompIS
 IS, target, dimension(1)          :: bcPointIS
 IS, pointer                       :: pBcCompIS(:)
 IS, pointer                       :: pBcPointIS(:)
 character(100)                    :: fname
 integer                           :: i, pst, pend
 integer                           :: vst, vend
 PetscInt                          :: size

 PetscInt                                :: nroots, nleaves
 PetscInt, target, dimension(1)          :: roots
 type(PetscSFNode), target, dimension(1) :: leaves
 PetscInt, pointer                       :: iroots(:)
 type(PetscSFNode), pointer              :: ileaves(:)



 call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 if(ierr.ne.0) then
   write(*,*) 'Fail to initialize PETSc ...'
   stop
 endif
 call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
 call MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
 write(*,*) 'MPI information ... ', myid, nproc


 fname = './square.msh'
 call DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, fname, PETSC_TRUE, dm_dumm, ierr);CHKERRA(ierr)

 !=========================================!
 ! define section                          !
 !=========================================!
 call DMGetDimension(dm_dumm, dim, ierr);CHKERRA(ierr)

 ! Create a scalar field, a vector field, and a surface vector field
 numComp(:) = 1 ! Note: dont know about exact meaning yet
 pNumComp => numComp
 do i=1,numFields*(dim+1)
    numDof(i) = 0
 end do

 ! define first field on vertices
 numDof(0*(dim+1)+1)     = 1 ! number of content in "vertex" field

 ! define second field on cells
 numDof(1*(dim+1)+dim+1) = 0 ! number of content in "cell" field

 ! define third field on faces
 numDof(2*(dim+1)+dim)   = 0 ! number of content in "face" field
 pNumDof => numDof

 ! Setup boundary conditions
 numBC = 1
 bcField(1) = 0
 pBcField => bcField
 call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr);CHKERRA(ierr)
 pBcCompIS => bcCompIS

 call DMGetStratumIS(dm_dumm, "depth", dim, bcPointIS(1),ierr);CHKERRA(ierr)
 pBcPointIS => bcPointIS

 ! Create a PetscSection with this data layout
 call DMSetNumFields(dm_dumm, numFields,ierr);CHKERRA(ierr)
 call DMPlexCreateSection(dm_dumm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr)
 CHKERRA(ierr)
 call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr)
 call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr)

 ! Name the Field variables
 call PetscSectionSetFieldName(section, 0, 'Node', ierr);CHKERRA(ierr)
 call DMSetLocalSection(dm_dumm, section, ierr);CHKERRA(ierr)
 call PetscSectionDestroy(section, ierr);CHKERRA(ierr)


 ! Tell PETSc to use natural indexing
 call DMSetUseNatural(dm_dumm, PETSC_TRUE, ierr);CHKERRA(ierr)

 ! Distribute mesh over processors
 if(nproc.gt.1) then
   call DMPlexDistribute(dm_dumm, 0, PETSC_NULL_SF, dm_geom, ierr);CHKERRA(ierr)
   call DMGetNaturalSF(dm_geom, sf_nat, ierr);CHKERRA(ierr)
   call PetscSFView(sf_nat, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)
 else
   call DMClone(dm_dumm, dm_geom, ierr);CHKERRA(ierr)
 endif

 ! extract natural order
 call PetscSFGetGraph(sf_nat, nroots, nleaves, iroots, ileaves, ierr);CHKERRA(ierr)
 ! check
 do i=1,nleaves
   write(*,*) myid, ileaves(i)%rank, ileaves(i)%index
 enddo

 call PetscFinalize(ierr)
 stop
end program test