BuildMap_SetInterface3D_t Subroutine

  subroutine BuildMap_SetInterface3D_t(this,globalMesh,setOfElem,setGlobalId, &
                                       setOffsetThisRank,nvarOfSet,N)
    !! Build the interface map from the global mesh and the set numbering, without
    !! communication beyond duplicating the communicator. Collective over
    !! globalMesh%decomp%mpiComm.
    !!
    !!  Input
    !!    - globalMesh : the decomposed global mesh
    !!    - setOfElem(1:nGlobalElem) : set id of every global element
    !!    - setGlobalId(1:nGlobalElem) : set-global id of every global element
    !!    - setOffsetThisRank(1:nSets) : number of set-k elements on lower ranks, so the
    !!      sub-mesh element of global element e is setGlobalId(e)-setOffsetThisRank(s(e))
    !!    - nvarOfSet(1:nSets) : number of variables in each set's trace
    !!    - N : polynomial degree; faces hold (N+1)*(N+1) nodes
    implicit none
    class(SetInterface3D_t),intent(inout) :: this
    type(Mesh3D),intent(in) :: globalMesh
    integer,intent(in) :: setOfElem(1:)
    integer,intent(in) :: setGlobalId(1:)
    integer,intent(in) :: setOffsetThisRank(1:)
    integer,intent(in) :: nvarOfSet(1:)
    integer,intent(in) :: N
    ! Local
    integer :: nSets,iEl,s1,eg,e2g,k,p,f,c,r,a,b,nn,ierror,np2,offset
    integer,allocatable :: counts(:),order(:),keys(:),faceOfSide(:,:)
    logical,allocatable :: isNbr(:)

    nSets = size(nvarOfSet)
    this%nSets = nSets
    this%N = N
    np2 = (N+1)*(N+1)
    this%nvarOfSet = nvarOfSet
    this%nvarMax = maxval(nvarOfSet)
    this%rankId = globalMesh%decomp%rankId
    this%mpiPrec = globalMesh%decomp%mpiPrec
    call mpi_comm_dup(globalMesh%decomp%mpiComm,this%comm,ierror)

    ! Count local interface faces per receiving set
    allocate(counts(1:nSets))
    counts = 0
    do iEl = 1,globalMesh%nElem
      eg = globalMesh%decomp%offsetElem(this%rankId+1)+iEl
      do s1 = 1,6
        e2g = globalMesh%sideInfo(3,s1,iEl)
        if(e2g <= 0) cycle
        if(setOfElem(e2g) == setOfElem(eg)) cycle
        counts(setOfElem(eg)) = counts(setOfElem(eg))+1
      enddo
    enddo

    allocate(this%setFaceStart(1:nSets+1))
    this%setFaceStart(1) = 1
    do k = 1,nSets
      this%setFaceStart(k+1) = this%setFaceStart(k)+counts(k)
    enddo
    this%nFaces = this%setFaceStart(nSets+1)-1
    nn = this%nFaces

    allocate(this%faceSet(1:nn),this%faceElem(1:nn),this%faceSide(1:nn), &
             this%faceFlip(1:nn),this%facePartnerSet(1:nn),this%facePartnerRank(1:nn), &
             this%facePartnerElem(1:nn),this%facePartnerSide(1:nn),this%faceSideId(1:nn), &
             this%facePartnerFace(1:nn))
    allocate(this%outbox(1:N+1,1:N+1,1:this%nvarMax,1:nn))
    allocate(this%ghost(1:N+1,1:N+1,1:this%nvarMax,1:nn))
    this%outbox = 0.0_prec
    this%ghost = 0.0_prec
    offset = globalMesh%decomp%offsetElem(this%rankId+1)
    allocate(faceOfSide(1:6,1:globalMesh%nElem))
    faceOfSide = 0

    ! Fill the faces, grouped by set, in global element order within each set
    counts = 0
    do iEl = 1,globalMesh%nElem
      eg = globalMesh%decomp%offsetElem(this%rankId+1)+iEl
      do s1 = 1,6
        e2g = globalMesh%sideInfo(3,s1,iEl)
        if(e2g <= 0) cycle
        k = setOfElem(eg)
        p = setOfElem(e2g)
        if(p == k) cycle
        f = this%setFaceStart(k)+counts(k)
        counts(k) = counts(k)+1
        faceOfSide(s1,iEl) = f
        this%faceSet(f) = k
        this%faceElem(f) = setGlobalId(eg)-setOffsetThisRank(k)
        this%faceSide(f) = s1
        this%facePartnerSide(f) = globalMesh%sideInfo(4,s1,iEl)/10
        this%faceFlip(f) = globalMesh%sideInfo(4,s1,iEl)-10*this%facePartnerSide(f)
        this%facePartnerSet(f) = p
        this%facePartnerRank(f) = globalMesh%decomp%elemToRank(e2g)
        this%faceSideId(f) = abs(globalMesh%sideInfo(2,s1,iEl))
        if(this%facePartnerRank(f) == this%rankId) then
          this%facePartnerElem(f) = setGlobalId(e2g)-setOffsetThisRank(p)
        else
          this%facePartnerElem(f) = 0
        endif
      enddo
    enddo
    deallocate(counts)

    ! Same-rank partner faces, through the partner's global element and side
    do f = 1,nn
      this%facePartnerFace(f) = 0
    enddo
    do iEl = 1,globalMesh%nElem
      do s1 = 1,6
        f = faceOfSide(s1,iEl)
        if(f == 0) cycle
        if(this%facePartnerRank(f) /= this%rankId) cycle
        e2g = globalMesh%sideInfo(3,s1,iEl)
        this%facePartnerFace(f) = faceOfSide(this%facePartnerSide(f),e2g-offset)
      enddo
    enddo
    deallocate(faceOfSide)

    ! Neighbour ranks, in increasing rank order
    allocate(isNbr(0:globalMesh%decomp%nRanks-1))
    isNbr = .false.
    do f = 1,nn
      if(this%facePartnerRank(f) /= this%rankId) isNbr(this%facePartnerRank(f)) = .true.
    enddo
    this%nNbr = count(isNbr)
    allocate(this%nbrRank(1:this%nNbr))
    c = 0
    do r = 0,globalMesh%decomp%nRanks-1
      if(isNbr(r)) then
        c = c+1
        this%nbrRank(c) = r
      endif
    enddo
    deallocate(isNbr)

    ! Faces per neighbour, ordered by global side id
    allocate(this%nbrFaceStart(1:this%nNbr+1))
    allocate(this%nbrFace(1:count(this%facePartnerRank /= this%rankId)))
    this%nbrFaceStart(1) = 1
    do c = 1,this%nNbr
      a = this%nbrFaceStart(c)
      b = a-1
      do f = 1,nn
        if(this%facePartnerRank(f) == this%nbrRank(c)) then
          b = b+1
          this%nbrFace(b) = f
        endif
      enddo
      this%nbrFaceStart(c+1) = b+1
      if(b > a) then
        allocate(order(1:b-a+1),keys(1:b-a+1))
        order = this%nbrFace(a:b)
        keys = this%faceSideId(order)
        call HeapSortByKey(keys,order)
        this%nbrFace(a:b) = order
        deallocate(order,keys)
      endif
    enddo

    ! Buffer offsets: a face sends nvar of its own set and receives nvar of the partner's
    allocate(this%sendStart(1:this%nNbr+1),this%recvStart(1:this%nNbr+1))
    this%sendStart(1) = 1
    this%recvStart(1) = 1
    do c = 1,this%nNbr
      this%sendStart(c+1) = this%sendStart(c)
      this%recvStart(c+1) = this%recvStart(c)
      do b = this%nbrFaceStart(c),this%nbrFaceStart(c+1)-1
        f = this%nbrFace(b)
        this%sendStart(c+1) = this%sendStart(c+1)+np2*nvarOfSet(this%faceSet(f))
        this%recvStart(c+1) = this%recvStart(c+1)+np2*nvarOfSet(this%facePartnerSet(f))
      enddo
    enddo
    allocate(this%sendBuf(1:this%sendStart(this%nNbr+1)-1))
    allocate(this%recvBuf(1:this%recvStart(this%nNbr+1)-1))

    ! Persistent requests: receives first, then sends
    allocate(this%requests(1:2*this%nNbr))
    do c = 1,this%nNbr
      call mpi_recv_init(this%recvBuf(this%recvStart(c)), &
                         this%recvStart(c+1)-this%recvStart(c),this%mpiPrec, &
                         this%nbrRank(c),SET_INTERFACE_TAG,this%comm, &
                         this%requests(c),ierror)
      call mpi_send_init(this%sendBuf(this%sendStart(c)), &
                         this%sendStart(c+1)-this%sendStart(c),this%mpiPrec, &
                         this%nbrRank(c),SET_INTERFACE_TAG,this%comm, &
                         this%requests(this%nNbr+c),ierror)
    enddo

  endsubroutine BuildMap_SetInterface3D_t