BuildHaloExchange_DomainDecomposition Subroutine

  subroutine BuildHaloExchange_DomainDecomposition(this,sideInfo,nElem,nSidesPerElem)
  !! Build the aggregated halo exchange tables for the mesh this domain
  !! decomposition belongs to.
  !!
  !! Enumerates every locally-owned side whose neighbor element resides on
  !! another rank, groups the sides by neighbor rank, and sorts each group by
  !! global side id. Both ranks sharing an interface enumerate the interface
  !! sides in the same order (the global side id is shared), so a packed send
  !! buffer on one rank and the packed receive buffer on the other line up
  !! without any per-message metadata. The (element,side,flip) triplets are
  !! copied to the device for the field pack/unpack kernels.
  !!
  !!  Input
  !!    - sideInfo : the mesh sideInfo array, (1:5, 1:nSidesPerElem, 1:nElem)
  !!    - nElem : number of local elements
  !!    - nSidesPerElem : 4 for 2-D (quadrilateral), 6 for 3-D (hexahedral)
    implicit none
    class(DomainDecomposition),intent(inout) :: this
    integer,intent(in) :: nElem
    integer,intent(in) :: nSidesPerElem
    integer,intent(in) :: sideInfo(1:5,1:nSidesPerElem,1:nElem)
    ! Local
    integer :: e1,s1,e2,r2,n,nhalo
    integer,allocatable :: elems(:),sides(:),flips(:),ranks(:)
    integer,allocatable,target :: halosides(:)
    integer(int64),allocatable :: keys(:)

    nhalo = 0
    do e1 = 1,nElem
      do s1 = 1,nSidesPerElem
        e2 = sideInfo(3,s1,e1)
        if(e2 > 0) then
          if(this%elemToRank(e2) /= this%rankId) then
            nhalo = nhalo+1
          endif
        endif
      enddo
    enddo

    this%halo_nsides = nhalo
    this%halo_nnbr = 0
    if(nhalo == 0) then
      this%halo_built = .true.
      return
    endif

    allocate(elems(1:nhalo),sides(1:nhalo),flips(1:nhalo), &
             ranks(1:nhalo),keys(1:nhalo))

    n = 0
    do e1 = 1,nElem
      do s1 = 1,nSidesPerElem
        e2 = sideInfo(3,s1,e1)
        if(e2 > 0) then
          r2 = this%elemToRank(e2)
          if(r2 /= this%rankId) then
            n = n+1
            elems(n) = e1
            sides(n) = s1
            flips(n) = sideInfo(4,s1,e1)-10*(sideInfo(4,s1,e1)/10)
            ranks(n) = r2
            ! Composite sort key: neighbor rank (major), global side id (minor).
            keys(n) = int(r2,int64)*2147483648_int64+ &
                      int(abs(sideInfo(2,s1,e1)),int64)
          endif
        endif
      enddo
    enddo

    call HaloKeySort(keys,elems,sides,flips,ranks,nhalo)

    ! Count the neighboring ranks and record the segment offsets
    this%halo_nnbr = 1
    do n = 2,nhalo
      if(ranks(n) /= ranks(n-1)) this%halo_nnbr = this%halo_nnbr+1
    enddo

    allocate(this%halo_rank(1:this%halo_nnbr))
    allocate(this%halo_offset(1:this%halo_nnbr+1))

    this%halo_nnbr = 1
    this%halo_rank(1) = ranks(1)
    this%halo_offset(1) = 0
    do n = 2,nhalo
      if(ranks(n) /= ranks(n-1)) then
        this%halo_nnbr = this%halo_nnbr+1
        this%halo_rank(this%halo_nnbr) = ranks(n)
        this%halo_offset(this%halo_nnbr) = n-1
      endif
    enddo
    this%halo_offset(this%halo_nnbr+1) = nhalo

    ! Flatten the (element,side,flip) triplets with 0-based element and side
    ! indices for the device pack/unpack kernels
    allocate(halosides(1:3*nhalo))
    do n = 1,nhalo
      halosides(3*n-2) = elems(n)-1
      halosides(3*n-1) = sides(n)-1
      halosides(3*n) = flips(n)
    enddo

    call gpuCheck(hipMalloc(this%halo_sides_gpu,sizeof(halosides)))
    call gpuCheck(hipMemcpy(this%halo_sides_gpu,c_loc(halosides), &
                            sizeof(halosides),hipMemcpyHostToDevice))

    deallocate(elems,sides,flips,ranks,keys,halosides)

    this%halo_built = .true.

  endsubroutine BuildHaloExchange_DomainDecomposition