GridSorting.f90

Go to the documentation of this file.

module gridsorting
  use kindmodule, only: i4b, dp, lgp
  use constantsmodule, only: dhalf, dzero
  use cellwithnbrsmodule, only: globalcelltype
  use mathutilmodule, only: is_close
  use basedismodule, only: dis_transform_xy
  implicit none
  private
 
  public :: quicksortgrid
 
contains
  ! Sort an array of integers
  subroutine quicksortgrid(array, arraySize, idxToGlobal, z_only)
    integer, intent(inout), dimension(:) :: array
    integer, intent(in) :: arraysize
    type(globalcelltype), dimension(:), pointer :: idxtoglobal
    logical(LGP) :: z_only !< only sort by z coordinate
    ! local
    integer :: qsort_threshold = 8
    include "qsort_inline.inc"
 
  contains
    subroutine init()
    end subroutine init
 
    ! Compare two grid cells, this doesn't work as
    ! smooth for staggered discretizations though...
    function lessthan(n, m) result(isLess)
      integer(I4B), intent(in) :: n
      integer(I4B), intent(in) :: m
      logical(LGP) :: isless
      ! local
      real(dp), dimension(3) :: xyz_n, xyz_m
 
      ! get coordinates as 3-vectors
      xyz_n = get_global_xyz(idxtoglobal(array(n)), z_only)
      xyz_m = get_global_xyz(idxtoglobal(array(m)), z_only)
 
      ! compare
      if (.not. is_close(xyz_n(3), xyz_m(3), 10 * epsilon(xyz_n(3)))) then
        isless = xyz_n(3) > xyz_m(3)
      else if (.not. is_close(xyz_n(2), xyz_m(2), 10 * epsilon(xyz_n(2)))) then
        isless = xyz_n(2) > xyz_m(2)
      else if (.not. is_close(xyz_n(1), xyz_m(1), 10 * epsilon(xyz_n(1)))) then
        isless = xyz_n(1) < xyz_m(1)
      else
        isless = .false.
      end if
 
    end function lessthan
 
    !> @brief Utility function to convert global cell
    !< id to global x,y,z coordinates
    function get_global_xyz(gc, use_only_z) result(global_xyz)
      type(globalcelltype) :: gc !< the global cell id
      logical(LGP) :: use_only_z !< only z coordinate is needed or available, skip transform
      real(dp), dimension(3) :: global_xyz !< return xyz
      ! local
      real(dp) :: x, y, z
      real(dp) :: xc, yc, xo, yo, angrot
 
      z = dhalf * (gc%v_model%dis_top%get(gc%index) + &
                   gc%v_model%dis_bot%get(gc%index))
 
      x = dzero
      y = dzero
      if (.not. use_only_z) then
        xc = gc%v_model%dis_xc%get(gc%index)
        yc = gc%v_model%dis_yc%get(gc%index)
        xo = gc%v_model%dis_xorigin%get()
        yo = gc%v_model%dis_yorigin%get()
        angrot = gc%v_model%dis_angrot%get()
        call dis_transform_xy(xc, yc, xo, yo, angrot, x, y)
      end if
 
      global_xyz = [x, y, z]
 
    end function get_global_xyz
 
    ! swap indices
    subroutine swap(a, b)
      integer, intent(in) :: a, b
      integer :: hold
 
      hold = array(a)
      array(a) = array(b)
      array(b) = hold
 
    end subroutine swap
 
    ! circular shift-right by one
    subroutine rshift(left, right)
      integer, intent(in) :: left, right
      integer :: hold
 
      hold = array(right)
      array(left + 1:right) = array(left:right - 1)
      array(left) = hold
 
    end subroutine rshift
  end subroutine quicksortgrid
end module gridsorting