program simpleTest2

  !use newsubs
  implicit none

#include <petsc-finclude/petsc.h90>

  ! Geometric dimension:
#define DIM 3

  MPI_Comm :: comm
  PetscMPIInt :: rank
  character(80) :: filename = "simpleblock-100.exo"
  PetscInt, parameter :: dof = 1
  PetscErrorCode :: ierr
  DM :: dm, dist_dm, ghost_dm
  DMLabel :: label
  PetscFV :: fvm
  PetscDS :: ds
  PetscInt :: num_ghost_cells
  PetscInt :: c, start_cell, end_cell, end_interior_cell
  PetscInt :: start_face, end_face, num_faces
  PetscInt :: overlap
  Vec :: X
  PetscScalar, pointer :: localx(:)
  DMLabel :: ghostLabel
  Vec :: cell_geom
  DM  :: geom_dm
  PetscSection :: geom_section, solution_section
  PetscInt :: cell_offset, geom_offset, ghost

  
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
  comm = PETSC_COMM_WORLD
  call MPI_Comm_rank(comm, rank, ierr); CHKERRQ(ierr)

  ! Set up mesh  
  call DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, dm, ierr); CHKERRQ(ierr)
  call DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE, ierr); CHKERRQ(ierr)
  call DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE, ierr); CHKERRQ(ierr)

  ! Distribute mesh
  overlap = 1
  call DMPlexDistribute(dm, overlap, PETSC_NULL_OBJECT, dist_dm, ierr)
  CHKERRQ(ierr)
  if (dist_dm /= 0) then
     call DMDestroy(dm, ierr); CHKERRQ(ierr)
     dm = dist_dm
  end if

  call DMSetFromOptions(dm, ierr); CHKERRQ(ierr)

! Create a label for the dm but don't assign it to any faces. This prevents 
! DMPlexConstructGhostCells from labelling the faces 
  call DMPlexCreateLabel(dm, "ClosedBoundaries", ierr); CHKERRQ(ierr)

  call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

! Construct ghost cells
! On a single processor there are no ghost cells but this labels all the outer
! boundaries as ghost faces which is good so that we don't try to work out the
! flux on those faces
! On two processeors this labels the two overlap cells as ghost cells and all
! the boundary faces as ghost faces.  
  
  call DMPlexConstructGhostCells(dm, "ClosedBoundaries", num_ghost_cells, ghost_dm, ierr)

  CHKERRQ(ierr)
  if (ghost_dm /= 0) then
     call DMDestroy(dm, ierr); CHKERRQ(ierr)
     dm = ghost_dm
  end if

  call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

  call DMPlexGetLabel(dm, "ghost", label, ierr); CHKERRQ(ierr)
  call PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE, ierr); CHKERRQ(ierr)
  call DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

! Set up discretization  
  call PetscFVCreate(comm, fvm, ierr); CHKERRQ(ierr)
  call PetscFVSetNumComponents(fvm, dof, ierr); CHKERRQ(ierr)
  call PetscFVSetSpatialDimension(fvm, DIM, ierr); CHKERRQ(ierr)
  call DMGetDS(dm, ds, ierr); CHKERRQ(ierr)
  call PetscDSAddDiscretization(ds, fvm, ierr); CHKERRQ(ierr)
  call PetscFVDestroy(fvm, ierr); CHKERRQ(ierr)

  ! Testing  
  call DMPlexGetHeightStratum(dm, 0, start_cell, end_cell, ierr); CHKERRQ(ierr)
  call DMPlexGetHybridBounds(dm, end_interior_cell, &
       PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr); CHKERRQ(ierr)
  call DMPlexGetHeightStratum(dm, 1, start_face, end_face, ierr); CHKERRQ(ierr)

! On two processors the end_cell = end_interior_cell values on both processors
! Is this correct? I guess we can identify the overlap cells as they have been labelled as ghost?  
  write(*,*) start_cell, end_cell, end_interior_cell, start_face, end_face

  call DMCreateGlobalVector(dm, X, ierr); CHKERRQ(ierr);
  call VecView(X, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

  call VecGetArrayF90(X, localx, ierr); CHKERRQ(ierr);
  call DMGetDefaultSection(dm, solution_section, ierr); CHKERRQ(ierr)

! We worked out its important not to write to the parts of arrays that correspond to overlap cells.
! Following ex11 it would be done something like this using C:
!
! ierr = DMPlexPointGlobalRef(dm,c,x,&xc);CHKERRQ(ierr);
! if (xc) {ierr = (*mod->solution)(mod,0.0,cg->centroid,xc,mod->solutionctx);CHKERRQ(ierr);}
!    
! And according to the documentation on DMPlexPointGlobalRef xc would be returned as null for
! the overlap blocks that don't belong to the processor
!
! We can achieve the same using the ghost labels on the overlap cells but is this a good
! idea? We'd prefer to follow the standard approach ie in ex11 but I don't think there's a
! clean way to do this in fortran as we can't pass pointers to structures...
! Could it be done with DMPlexPointGlobalRef? If so, any ideas?
    
  do c = start_cell, end_interior_cell-1    
! ???  call DMPlexPointGlobalRef(dm, c, localx, someref)???
    call DMLabelGetValue(label, c, ghost, ierr); CHKERRQ(ierr)
    if (ghost < 0) then
       call PetscSectionGetOffset(solution_section, c, cell_offset, ierr); CHKERRQ(ierr)
       cell_offset = cell_offset + 1
       localx(cell_offset) = 100.0
    end if
 end do
   
  call VecRestoreArrayF90(X, localx, ierr); CHKERRQ(ierr);

  call VecView(X, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)
  call VecDestroy(X, ierr); CHKERRQ(ierr)

  call DMDestroy(dm, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program simpleTest2