BuildSetNumbering Subroutine

  subroutine BuildSetNumbering(globalMesh,setOfElem,nSets,setGlobalId,setOffsets)
    !! Number the elements of every set and partition each set over the ranks.
    !!
    !! For each global element e with k = s(e), setGlobalId(e) = g_k(e), its id within
    !! set k (see the module documentation). For each set k, setOffsets(1:nRanks+1,k)
    !! are the decomposition offsets of the set's sub-mesh: rank r owns set-global ids
    !! setOffsets(r+1,k)+1 .. setOffsets(r+2,k), which are exactly the set-k elements in
    !! rank r's range of the global mesh. A rank that owns no element of set k has
    !! setOffsets(r+1,k) = setOffsets(r+2,k). Rank-local; every rank gets the same result.
    !!
    !!  Input
    !!    - globalMesh : the decomposed global mesh
    !!    - setOfElem(1:nGlobalElem) : set id of each global element, validated
    !!    - nSets : number of sets
    !!  Output
    !!    - setGlobalId(1:nGlobalElem) : set-global id of each global element
    !!    - setOffsets(1:nRanks+1,1:nSets) : per-set decomposition offsets
    implicit none
    type(Mesh3D),intent(in) :: globalMesh
    integer,intent(in) :: setOfElem(1:)
    integer,intent(in) :: nSets
    integer,intent(out) :: setGlobalId(1:)
    integer,intent(out) :: setOffsets(1:,1:)
    ! Local
    integer :: e,k,r,nRanks
    integer,allocatable :: nInSet(:)

    nRanks = globalMesh%decomp%nRanks
    allocate(nInSet(1:nSets))
    nInSet = 0
    setOffsets(1:nRanks+1,1:nSets) = 0

    do r = 1,nRanks
      do e = globalMesh%decomp%offsetElem(r)+1,globalMesh%decomp%offsetElem(r+1)
        k = setOfElem(e)
        nInSet(k) = nInSet(k)+1
        setGlobalId(e) = nInSet(k)
      enddo
      setOffsets(r+1,1:nSets) = nInSet(1:nSets)
    enddo

    deallocate(nInSet)

  endsubroutine BuildSetNumbering