module PETSc_Tools

#include "petsc/finclude/petscdmplex.h"
#include "petsc/finclude/petscdmlabel.h"

 use petscsys
 use petscdmplex
 use IO_Treat, only: Outputs

 implicit none

 public ! default as public

 !----------------------------------!
 ! General PETSc variables          !
 !----------------------------------!
 PetscErrorCode      :: ierr
 PetscMPIInt         :: myid, nproc, icpu
 PetscInt            :: size

 !----------------------------------!
 ! PETSc section and output field   !
 !----------------------------------!
 DMLabel, pointer                  :: nolabel(:) => NULL()
 Vec                               :: node          ! node looping vector
 PetscSection                      :: section
 PetscInt                          :: numFields = 3 ! Note: numFields does not mean "number of variables"
                                                    ! Example: 
                                                    ! numFields=1: put data on vertices
                                                    ! numFields=2: put data on cells
                                                    ! numFields=3: put data on faces
 PetscInt                          :: numBC     = 1
 PetscInt, target, dimension(3)    :: numComp ! dont fix size. Its okay for Finite-Volume
 PetscInt, target, dimension(12)   :: numDof  ! dont fix size. Its okay for Finite-Volume
 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(:)
 PetscViewer                       :: viewer

 !----------------------------------!
 ! DMPlex variables to handle grid  !
 !----------------------------------!
 DM                  :: dm_geom               ! dmplex object encapsulates loaded mesh
 Vec                 :: coord                 ! coordinate array from PETSc
 PetscInt            :: depth, dim            ! depth and dimension in spatial grid
 PetscInt            :: pst, pend             ! start & end index for entire topological points
 PetscInt            :: vst, vend             ! start & end index for verticies
 PetscInt            :: est, eend             ! start & end index for edges
 PetscInt            :: cst, cend             ! start & end index for cells 
 PetscInt            :: fst, fend             ! start & end index for faces 
 PetscInt            :: nnode                 ! number of vertices 
 PetscInt            :: nedge                 ! number of edges
 PetscInt            :: ncell                 ! number of cells
! PetscInt            :: nface                 ! number of faces
 PetscInt            :: nebnd                 ! number of boundary edge (physical)
 PetscInt            :: nebnda                ! number of boundary edge (physical & parallel)
 PetscInt            :: nebndp                ! number of boundary edge (parallel)
 PetscInt            :: neinn                 ! number of inner edge (physical & parallel boundaries excluded)
 DMLabel             :: label_physical        ! label for physical boundary (exclude parallel boundary)
 DMLabel             :: label_para            ! label for parallel boundary 
 DMLabel             :: label_all             ! label for all kinds of boundary (both physical and parallel)
 DMLabel             :: label_cell            ! label for cell shape
 PetscInt, parameter :: marker_boundary = 100 ! marker that includes both physical and parallel boundaries
 PetscInt, parameter :: marker_para     = 200 ! marker for parallel boundary

 !---------------------------------------!
 ! PETSc objects for halo communication  !
 !---------------------------------------!
 PetscSF             :: sf!, point_sf, section_sf
 PetscSection        :: s
 type(tVec)          :: vloc ! local vector
 type(tVec)          :: vglo ! global vector

 !------------------------------------------!
 ! Dummy variables for vector manipulation  !
 !------------------------------------------!
  !-----------------------!
  ! general for 2D and 3D !
  !-----------------------!
 PetscReal, target, dimension(1) :: sol_tmp
 PetscReal, pointer              :: pSol(:) !=> sol_tmp
 PetscReal, target, dimension(1) :: xx, yy, zz
 PetscReal, pointer              :: pxx(:)  !=> xx
 PetscReal, pointer              :: pyy(:)  !=> yy
 PetscReal, pointer              :: pzz(:)  !=> zz
 double precision, allocatable   :: sol(:)  ! temporary solution vector. 

  !-----------------------!
  ! for 2D                !
  !-----------------------!
 PetscInt, target, dimension(4)  :: cone_2D          ! number 4 covers triangle and quadrilateral cells in 2D
 PetscInt, pointer               :: pCone_2D(:)      !=> cone_2D
 PetscInt, pointer               :: pCone_tmp_2D(:)  !=> cone_2D
 PetscInt, target, dimension(10) :: supp_2D          ! number 10 denotes upper limit of supporting edges at given vertex in 2D
 PetscInt, pointer               :: pSupp_2D(:)      !=> supp_2D
 PetscReal, target, dimension(2) :: cent_2D, norm_2D ! number 2 means x and y coordinates
 PetscReal, pointer              :: pCent_2D(:)      !=> cent_2D
 PetscReal, pointer              :: pNorm_2D(:)      !=> norm_2D

contains

 !----------------------------------!
 subroutine PETSc_Output_Field()
  ! print out field data into a vtu file

  use Utils_Parameters

  implicit none

  integer          :: i
  type(Outputs)    :: Output
  character(str_l) :: fname


  ! 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
                              ! if this number is equal to 1, scalar will be printed
                              ! if this number is lager than 1, array will be printed

  ! define second field on cells
  numDof(1*(dim+1)+dim+1) = 0 ! number of content in "cell" field
                              ! Note: even without defining cell field, 
                              !       Rank field still printed out

  ! define third field on faces
  numDof(2*(dim+1)+dim)   = 0 ! number of content in "face" field
                              ! Note: face field will not shown in VTK
  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

  ! Note : should turn on marker label using option (see run file)
  !        using other labels failed
  call DMGetStratumIS(dm_geom, "marker", 1, bcPointIS(1),ierr);CHKERRA(ierr)

  pBcPointIS => bcPointIS

  ! Create a PetscSection with this data layout
  ! nolabel(:) => NULL()
  call DMSetNumFields(dm_geom, numFields,ierr);CHKERRA(ierr)
  call DMPlexCreateSection(dm_geom,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) ! value at node
  call DMSetLocalSection(dm_geom, section, ierr);CHKERRA(ierr)
  call PetscSectionDestroy(section, ierr);CHKERRA(ierr)

  ! Create a Vec with this layout of section
  ! Note : vector 'node' will include entire fields defined above. 
  call DMGetLocalVector(dm_geom, node, ierr);CHKERRA(ierr) 

  ! Note : Local vector size depends on the way section defined
  call VecGetSize(node, size, ierr);CHKERRA(ierr)


  do i=0,size-1 ! vertex first and then face or vertex : according to your section definition
    ! example - 1 : fill vertex vector and face factor separately
    ! if(i.le.nnode-1) then ! fill node
    !   ! do whatever you want, and then you fill solution vector as below before print it out
    !   ! using section layout we defined above
    !   call VecSetValues(node, 1, i, one, INSERT_VALUES, ierr);CHKERRA(ierr)
    ! else ! fill face (or cell center value)
    !   call VecSetValues(node, 1, i, two, INSERT_VALUES, ierr);CHKERRA(ierr)
    ! endif

    ! example - 2 : pass arbitrary non-petsc solution vector to petsc data layout
    pSol = sol(i+1) 
    call VecSetValues(node, 1, i, pSol, INSERT_VALUES, ierr);CHKERRA(ierr)
  enddo

  ! create viewer
  call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr)
 
  ! set viewer type
  call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr)
 
  ! set output file name for viewer
  call Output%Path
  fname = Output%Output_Path_PETSc_Field
  call PetscViewerFileSetName(viewer, fname, ierr);CHKERRA(ierr) 
  !call PetscViewerFileSetName(viewer, 'sol.vtu', ierr);CHKERRA(ierr) 
 
  ! print dm object to the file
  call VecView(node, viewer, ierr);CHKERRA(ierr)
 
  ! destroy
  call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr)
  call DMRestoreLocalVector(dm_geom, node, ierr);CHKERRA(ierr)
  call DMDestroy(dm_geom, ierr);CHKERRA(ierr)

 end subroutine PETSc_Output_Field

end module PETSc_Tools
